diff --git a/Helicity/Vertex/Scalar/FFSVertex.cc b/Helicity/Vertex/Scalar/FFSVertex.cc --- a/Helicity/Vertex/Scalar/FFSVertex.cc +++ b/Helicity/Vertex/Scalar/FFSVertex.cc @@ -1,129 +1,121 @@ // -*- C++ -*- // // FFSVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FFSVertex class. // #include "FFSVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace ThePEG::Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGFFSVertex("ThePEG::FFSVertex", "libThePEG.so"); void FFSVertex::Init() { static ClassDocumentation documentation ("The FFSVertex class is the implementation of the FFS" "vertex. All such vertices shoud inherit from it"); } // evaluate the full vertex Complex FFSVertex::evaluate(Energy2 q2, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, const ScalarWaveFunction & sca) { // calculate the couplings setCoupling(q2,sp.particle(),sbar.particle(),sca.particle()); Complex vertex( _left*(sbar.s1()*sp.s1()+sbar.s2()*sp.s2()) +_right*(sbar.s3()*sp.s3()+sbar.s4()*sp.s4()) ); // final factors return Complex(0.,1.)*norm()*sca.wave()*vertex; } // off-shell scalar ScalarWaveFunction FFSVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, complex mass, complex width) { // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+sp.momentum(); // first calculate the couplings setCoupling(q2,sp.particle(),sbar.particle(),out); Energy2 p2 = pout.m2(); Complex fact = -norm()*propagator(iopt,p2,out,mass,width); Complex output = _left*(sbar.s1()*sp.s1()+sbar.s2()*sp.s2()) +_right*(sbar.s3()*sp.s3()+sbar.s4()*sp.s4()); // final factors and output output*=fact; return ScalarWaveFunction(pout,out,output); } // off-shell spinor SpinorWaveFunction FFSVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out, const SpinorWaveFunction & sp, const ScalarWaveFunction & sca, complex mass, complex width) { // work out the momentum of the off-shell particle Lorentz5Momentum pout = sp.momentum()+sca.momentum(); // first calculate the couplings setCoupling(q2,sp.particle(),out,sca.particle()); Energy2 p2 = pout.m2(); Complex fact = -norm()*sca.wave()*propagator(iopt,p2,out,mass,width); Complex ii(0.,1.); // useful combinations of the momenta if(mass.real() < ZERO) mass = out->mass(); complex p1p2 = pout.x()+ii*pout.y(); complex p1m2 = pout.x()-ii*pout.y(); Complex s1(0.),s2(0.),s3(0.),s4(0.); LorentzSpinor spt = sp.wave(); complex p0p3=pout.e()+pout.z(); complex p0m3=pout.e()-pout.z(); - s1 = UnitRemoval::InvE * - fact*( _left*mass*spt.s1()+_right*(p0m3*spt.s3()-p1m2*spt.s4())); - s2 = UnitRemoval::InvE * - fact*( _left*mass*spt.s2()+_right*(p0p3*spt.s4()-p1p2*spt.s3())); - s3 = UnitRemoval::InvE * - fact*(_right*mass*spt.s3()+ _left*(p0p3*spt.s1()+p1m2*spt.s2())); - s4 = UnitRemoval::InvE * - fact*(_right*mass*spt.s4()+ _left*(p0m3*spt.s2()+p1p2*spt.s1())); + s1 = Complex(UnitRemoval::InvE * fact*( _left*mass*spt.s1()+_right*(p0m3*spt.s3()-p1m2*spt.s4()))); + s2 = Complex(UnitRemoval::InvE * fact*( _left*mass*spt.s2()+_right*(p0p3*spt.s4()-p1p2*spt.s3()))); + s3 = Complex(UnitRemoval::InvE * fact*(_right*mass*spt.s3()+ _left*(p0p3*spt.s1()+p1m2*spt.s2()))); + s4 = Complex(UnitRemoval::InvE * fact*(_right*mass*spt.s4()+ _left*(p0m3*spt.s2()+p1p2*spt.s1()))); return SpinorWaveFunction(pout,out,s1,s2,s3,s4); } // off-shell SpinorBar SpinorBarWaveFunction FFSVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const SpinorBarWaveFunction & sbar, const ScalarWaveFunction & sca, complex mass, complex width) { // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+sca.momentum(); // first calculate the couplings setCoupling(q2,out,sbar.particle(),sca.particle()); Energy2 p2 = pout.m2(); Complex fact = -norm()*sca.wave()*propagator(iopt,p2,out,mass,width); Complex ii(0.,1.); // momentum components if(mass.real() < ZERO) mass = out->mass(); complex p1p2 = pout.x()+ii*pout.y(); complex p1m2 = pout.x()-ii*pout.y(); // complex numbers for the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); LorentzSpinorBar sbart=sbar.wave(); complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); - s1 = UnitRemoval::InvE * - fact*( mass*_left*sbart.s1()-_right*(p0p3*sbart.s3()+p1p2*sbart.s4())); - s2 = UnitRemoval::InvE * - fact*( mass*_left*sbart.s2()-_right*(p1m2*sbart.s3()+p0m3*sbart.s4())); - s3 = UnitRemoval::InvE * - fact*(mass*_right*sbart.s3()- _left*(p0m3*sbart.s1()-p1p2*sbart.s2())); - s4 = UnitRemoval::InvE * - fact*(mass*_right*sbart.s4()+ _left*(p1m2*sbart.s1()-p0p3*sbart.s2())); + s1 = Complex(UnitRemoval::InvE * fact*( mass*_left*sbart.s1()-_right*(p0p3*sbart.s3()+p1p2*sbart.s4()))); + s2 = Complex(UnitRemoval::InvE * fact*( mass*_left*sbart.s2()-_right*(p1m2*sbart.s3()+p0m3*sbart.s4()))); + s3 = Complex(UnitRemoval::InvE * fact*(mass*_right*sbart.s3()- _left*(p0m3*sbart.s1()-p1p2*sbart.s2()))); + s4 = Complex(UnitRemoval::InvE * fact*(mass*_right*sbart.s4()+ _left*(p1m2*sbart.s1()-p0p3*sbart.s2()))); return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4); } diff --git a/Helicity/Vertex/Scalar/GeneralVSSVertex.cc b/Helicity/Vertex/Scalar/GeneralVSSVertex.cc --- a/Helicity/Vertex/Scalar/GeneralVSSVertex.cc +++ b/Helicity/Vertex/Scalar/GeneralVSSVertex.cc @@ -1,84 +1,84 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the GeneralVSSVertex class. // #include "GeneralVSSVertex.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeHelicityGeneralVSSVertex("ThePEG::GeneralVSSVertex", "libThePEG.so"); void GeneralVSSVertex::Init() { static ClassDocumentation documentation ("The GeneralVSSVertex class implements a general form of the vector-scalar-scalar interaction"); } // evaluate the vertex Complex GeneralVSSVertex::evaluate(Energy2 q2, const VectorWaveFunction & vec, const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2) { // calculate the coupling setCoupling(q2,vec.particle(),sca1.particle(),sca2.particle()); // calculate the vertex return UnitRemoval::InvE * -Complex(0.,1.) * norm() * sca1.wave()*sca2.wave()* vec.wave().dot(a_*sca1.momentum()+b_*sca2.momentum()); } // off-shell vector VectorWaveFunction GeneralVSSVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2, complex mass, complex width) { // outgoing momentum Lorentz5Momentum pout(sca1.momentum()+sca2.momentum()); // calculate the coupling setCoupling(q2,out,sca1.particle(),sca2.particle()); // mass and width if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); // calculate the prefactor Energy2 p2 = pout.m2(); Complex fact = norm()*sca1.wave()*sca2.wave()*propagator(iopt,p2,out,mass,width); // compute the vector LorentzPolarizationVector vec = -UnitRemoval::InvE * fact * (b_*sca2.momentum()+a_*sca1.momentum()); // massive outgoing vector if(mass.real()!=ZERO) { // first the dot product for the second term complex dot = vec*pout/mass2; vec -= dot*pout; } return VectorWaveFunction(pout,out,vec); } // return an off-shell scalar ScalarWaveFunction GeneralVSSVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const VectorWaveFunction & vec, const ScalarWaveFunction & sca, complex mass, complex width ) { // momentum of the particle Lorentz5Momentum pout = sca.momentum()+vec.momentum(); // calculate the coupling setCoupling(q2,vec.particle(),sca.particle(),out); // calculate the prefactor Energy2 p2 = pout.m2(); Complex fact = norm()*sca.wave()*propagator(iopt,p2,out,mass,width); // compute the wavefunction - fact = UnitRemoval::InvE * fact*vec.wave().dot(a_*sca.momentum()-b_*pout); + fact = Complex(UnitRemoval::InvE * fact*vec.wave().dot(a_*sca.momentum()-b_*pout)); return ScalarWaveFunction(pout,out,fact); } diff --git a/Helicity/Vertex/Scalar/GeneralVVSVertex.cc b/Helicity/Vertex/Scalar/GeneralVVSVertex.cc --- a/Helicity/Vertex/Scalar/GeneralVVSVertex.cc +++ b/Helicity/Vertex/Scalar/GeneralVVSVertex.cc @@ -1,125 +1,125 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the GeneralVVSVertex class. // #include "GeneralVVSVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Helicity/epsilon.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGGeneralVVSVertex("ThePEG::GeneralVVSVertex", "libThePEG.so"); void GeneralVVSVertex::Init() { static ClassDocumentation documentation ("The GeneralVVSVertex class implements a general form of" " the vector-vector-scalar interaction"); } Complex GeneralVVSVertex::evaluate(Energy2 q2,const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const ScalarWaveFunction & sca) { Lorentz5Momentum pSca = sca.momentum(); Lorentz5Momentum pvec1 = vec1.momentum(); Lorentz5Momentum pvec2 = vec2.momentum(); // calculate kinematics if(kinematics()) calculateKinematics(pSca,pvec1,pvec2); // calculate coupling setCoupling(q2, vec1.particle(), vec2.particle(), sca.particle()); Complex e1e2(vec1.wave().dot(vec2.wave())); complex e1p1(vec1.wave().dot(pvec1)); complex e1p2(vec1.wave().dot(pvec2)); complex e2p1(vec2.wave().dot(pvec1)); complex e2p2(vec2.wave().dot(pvec2)); complex p1p2(invariant(1,2)); LorentzPolarizationVectorE eps = epsilon(vec1.wave(),vec2.wave(),pvec2); complex p1Ep2 = eps.dot(pvec1); Complex output = UnitRemoval::InvE2 * (_a00*e1e2*p1p2 + _aEp*p1Ep2 + _a11*e1p1*e2p1 + _a12*e1p1*e2p2 + _a21*e1p2*e2p1+ _a22*e1p2*e2p2); return -norm()*Complex(0.,1.) * sca.wave() * output; } ScalarWaveFunction GeneralVVSVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, complex mass, complex width) { // pointers to the particle data objects tcPDPtr Pvec1(vec1.particle()); tcPDPtr Pvec2(vec2.particle()); Lorentz5Momentum pvec1 = vec1.momentum(); Lorentz5Momentum pvec2 = vec2.momentum(); Lorentz5Momentum pout = pvec1 + pvec2; pout.rescaleMass(); // calculate kinematics if needed if(kinematics()) calculateKinematics(-pout,pvec1,pvec2); // calculate coupling setCoupling(q2,Pvec1,Pvec2,out); // propagator Complex prop = propagator(iopt,pout.m2(),out,mass,width); // lorentz part Complex e1e2(vec1.wave().dot(vec2.wave())); complex e1p1(vec1.wave().dot(pvec1)); complex e1p2(vec1.wave().dot(pvec2)); complex e2p1(vec2.wave().dot(pvec1)); complex e2p2(vec2.wave().dot(pvec2)); complex p1p2(invariant(1,2)); LorentzPolarizationVectorE eps = epsilon(vec1.wave(),vec2.wave(),pvec2); complex p1Ep2 = eps.dot(pvec1); Complex output = UnitRemoval::InvE2 * (_a00*e1e2*p1p2 + _aEp*p1Ep2 + _a11*e1p1*e2p1 + _a12*e1p1*e2p2 + _a21*e1p2*e2p1+ _a22*e1p2*e2p2); output *= norm()*prop; return ScalarWaveFunction(pout,out,output); } VectorWaveFunction GeneralVVSVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const VectorWaveFunction & vec, const ScalarWaveFunction & sca, complex mass, complex width) { Lorentz5Momentum pSca = sca.momentum(); Lorentz5Momentum pvec1 = vec.momentum()+sca.momentum(); Lorentz5Momentum pvec2 = vec.momentum(); // calculate kinematics if(kinematics()) calculateKinematics(pSca,-pvec1,pvec2); // calculate coupling setCoupling(q2, out, vec.particle(), sca.particle()); // prefactor Energy2 p2 = pvec1.m2(); if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); Complex fact = -norm()* sca.wave() * propagator(iopt,p2,out,mass,width); // vertex as polarization vector complex e2p1(-vec.wave().dot(pvec1)); complex e2p2(vec.wave().dot(pvec2)); complex p1p2(invariant(1,2)); - LorentzPolarizationVector pv = (UnitRemoval::InvE2*_a00*p1p2*vec.wave() - - UnitRemoval::InvE2*_a11*e2p1*pvec1 - - UnitRemoval::InvE2*_a12*e2p2*pvec1 + - UnitRemoval::InvE2*_a21*e2p1*pvec2 + - UnitRemoval::InvE2*_a22*e2p2*pvec2 + - UnitRemoval::InvE2*_aEp*epsilon(pvec1,vec.wave(),pvec2)); + LorentzPolarizationVector pv = (LorentzPolarizationVector(UnitRemoval::InvE2*_a00*p1p2*vec.wave()) - + LorentzPolarizationVector(UnitRemoval::InvE2*_a11*e2p1*pvec1) - + LorentzPolarizationVector(UnitRemoval::InvE2*_a12*e2p2*pvec1) + + LorentzPolarizationVector(UnitRemoval::InvE2*_a21*e2p1*pvec2) + + LorentzPolarizationVector(UnitRemoval::InvE2*_a22*e2p2*pvec2) + + LorentzPolarizationVector(UnitRemoval::InvE2*_aEp*epsilon(pvec1,vec.wave(),pvec2))); // evaluate the wavefunction LorentzPolarizationVector vect; // massless case if(mass.real()==ZERO) { vect = fact*pv; } // massive case else { complex dot = pv.dot(pvec1)/mass2; vect = fact*(pv-dot*pvec1); } return VectorWaveFunction(pvec1,out,vect); } diff --git a/Helicity/Vertex/Scalar/VVVSVertex.cc b/Helicity/Vertex/Scalar/VVVSVertex.cc --- a/Helicity/Vertex/Scalar/VVVSVertex.cc +++ b/Helicity/Vertex/Scalar/VVVSVertex.cc @@ -1,197 +1,197 @@ // -*- C++ -*- // // VVVSVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the VVVSVertex class. // #include "VVVSVertex.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/epsilon.h" using namespace ThePEG; using namespace Helicity; void VVVSVertex::persistentOutput(PersistentOStream & os) const { os << scalar_; } void VVVSVertex::persistentInput(PersistentIStream & is, int) { is >> scalar_; } DescribeAbstractClass describeAbstractVVVSVertex("Helicity::VVVSVertex","libThePEG.so"); void VVVSVertex::Init() { static ClassDocumentation documentation ("The VVVSVertex class implements the helicity amplitude" "calculations for the triple gauge boson scalar vertex. Any " "implementation of such a vertex should inherit from in and implement" " the virtual setCoupling member to calculate the coupling." "As this vertex is not present in renormalisabe theories we assume the form from Higgs" "effective theory, i.e. the triple vector vertex multiplied by the scale wavefunction."); } // evaluate the vertex Complex VVVSVertex::evaluate(Energy2 q2, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const VectorWaveFunction & vec3, const ScalarWaveFunction & sca) { // calculate the coupling setCoupling(q2,vec1.particle(),vec2.particle(),vec3.particle(),sca.particle()); // the scalar form ofr the vertex if(scalar_) { complex alpha1(ZERO); // decide if we need to use special treatment to avoid gauge cancelations // first vector if(abs(vec1.t())!=0.) { if(abs(vec1.t())>0.1*max( max(abs(vec1.x()),abs(vec1.y())),abs(vec1.z()))) alpha1=vec1.e()/vec1.t(); } // second vector if(abs(vec2.t())!=0.) { if(abs(vec2.t())>0.1*max( max(abs(vec2.x()),abs(vec2.y())),abs(vec2.z()))) alpha1=vec2.e()/vec2.t(); } // third vector if(abs(vec3.t())!=0.) { if(abs(vec3.t())>0.1*max( max(abs(vec3.x()),abs(vec3.y())),abs(vec3.z()))) alpha1=vec3.e()/vec3.t(); } // dot products of the polarization vectors Complex dot12 = vec1.wave().dot(vec2.wave()); Complex dot13 = vec1.wave().dot(vec3.wave()); Complex dot23 = vec3.wave().dot(vec2.wave()); // dot products of polarization vectors and momentum complex dotp13 = vec3.wave().dot(LorentzPolarizationVectorE(vec1.momentum()) - alpha1 * vec1.wave()); complex dotp23 = vec3.wave().dot(LorentzPolarizationVectorE(vec2.momentum()) - alpha1 * vec2.wave()); complex dotp21 = vec1.wave().dot(LorentzPolarizationVectorE(vec2.momentum()) - alpha1 * vec2.wave()); complex dotp31 = vec1.wave().dot(LorentzPolarizationVectorE(vec3.momentum()) - alpha1 * vec3.wave()); complex dotp32 = vec2.wave().dot(LorentzPolarizationVectorE(vec3.momentum()) - alpha1 * vec3.wave()); complex dotp12 = vec2.wave().dot(LorentzPolarizationVectorE(vec1.momentum()) - alpha1 * vec1.wave()); // finally calculate the vertex return Complex(0.,1.)*norm()*UnitRemoval::InvE*sca.wave()* (dot12*(dotp13-dotp23)+dot23*(dotp21-dotp31)+dot13*(dotp32-dotp12)); } else { LorentzPolarizationVector eps = epsilon(vec1.wave(),vec2.wave(),vec3.wave()); complex dot = eps.dot(vec1.momentum()+vec2.momentum()+vec3.momentum()); return norm()*Complex(0.,1.) * sca.wave() * dot *UnitRemoval::InvE; } } // off-shell vector VectorWaveFunction VVVSVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const ScalarWaveFunction & sca, complex mass, complex width) { // output momenta Lorentz5Momentum pout =vec1.momentum()+vec2.momentum()+sca.momentum(); // calculate the coupling setCoupling(q2,out,vec1.particle(),vec2.particle(),sca.particle()); // prefactor Energy2 p2 = pout.m2(); Complex fact = norm()*propagator(iopt,p2,out,mass,width); if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); LorentzPolarizationVector vect; if(scalar_) { // dot products we need Complex dot12 = vec1.wave().dot(vec2.wave()); complex dota = vec1.wave().dot(pout+vec2.momentum()); complex dotb = vec2.wave().dot(pout+vec1.momentum()); // compute the polarization vector vect = UnitRemoval::InvE*fact*sca.wave()* (dot12*(vec1.momentum()-vec2.momentum())-dotb*vec1.wave()+dota*vec2.wave()); } else { vect = -UnitRemoval::InvE*fact*sca.wave()* epsilon(vec1.momentum()+vec2.momentum()+pout,vec1.wave(),vec2.wave()); } // scalar piece for massive case if(mass.real()!=ZERO) { complex dot = vect.dot(pout)/mass2; vect -= dot*pout; } return VectorWaveFunction(pout,out,vect); } // evaluate the off-shell scalar wavefunction ScalarWaveFunction VVVSVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const VectorWaveFunction & vec3, complex mass, complex width) { // calculate the coupling setCoupling(q2,vec1.particle(),vec2.particle(),vec3.particle(),out); // outgoing momentum Lorentz5Momentum pout = vec1.momentum()+vec2.momentum()+vec3.momentum(); Complex fact(0.); // wavefunction Energy2 p2 = pout.m2(); if(scalar_) { // calculate the triple vector piece complex alpha1(ZERO); // decide if we need to use special treatment to avoid gauge cancelations // first vector if(abs(vec1.t())!=0.) { if(abs(vec1.t())>0.1*max( max(abs(vec1.x()),abs(vec1.y())),abs(vec1.z()))) alpha1=vec1.e()/vec1.t(); } // second vector if(abs(vec2.t())!=0.) { if(abs(vec2.t())>0.1*max( max(abs(vec2.x()),abs(vec2.y())),abs(vec2.z()))) alpha1=vec2.e()/vec2.t(); } // third vector if(abs(vec3.t())!=0.) { if(abs(vec3.t())>0.1*max( max(abs(vec3.x()),abs(vec3.y())),abs(vec3.z()))) alpha1=vec3.e()/vec3.t(); } // dot products of the polarization vectors Complex dot12 = vec1.wave().dot(vec2.wave()); Complex dot13 = vec1.wave().dot(vec3.wave()); Complex dot23 = vec3.wave().dot(vec2.wave()); // dot products of polarization vectors and momentum complex dotp13 = vec3.wave().dot(LorentzPolarizationVectorE(vec1.momentum()) - alpha1 * vec1.wave()); complex dotp23 = vec3.wave().dot(LorentzPolarizationVectorE(vec2.momentum()) - alpha1 * vec2.wave()); complex dotp21 = vec1.wave().dot(LorentzPolarizationVectorE(vec2.momentum()) - alpha1 * vec2.wave()); complex dotp31 = vec1.wave().dot(LorentzPolarizationVectorE(vec3.momentum()) - alpha1 * vec3.wave()); complex dotp32 = vec2.wave().dot(LorentzPolarizationVectorE(vec3.momentum()) - alpha1 * vec3.wave()); complex dotp12 = vec2.wave().dot(LorentzPolarizationVectorE(vec1.momentum()) - alpha1 * vec1.wave()); // finally calculate and return the wavefunction - fact = -norm()*UnitRemoval::InvE*propagator(iopt,p2,out,mass,width)* - (dot12*(dotp13-dotp23)+dot23*(dotp21-dotp31)+dot13*(dotp32-dotp12)); + fact = -Complex(norm()*UnitRemoval::InvE*propagator(iopt,p2,out,mass,width)* + (dot12*(dotp13-dotp23)+dot23*(dotp21-dotp31)+dot13*(dotp32-dotp12))); } else { LorentzPolarizationVector eps = epsilon(vec1.wave(),vec2.wave(),vec3.wave()); complex dot = eps.dot(pout); - fact = norm()*UnitRemoval::InvE*propagator(iopt,p2,out,mass,width)*dot; + fact = Complex(norm()*UnitRemoval::InvE*propagator(iopt,p2,out,mass,width)*dot); } return ScalarWaveFunction(pout,out,fact); } diff --git a/Helicity/Vertex/Tensor/FFTVertex.cc b/Helicity/Vertex/Tensor/FFTVertex.cc --- a/Helicity/Vertex/Tensor/FFTVertex.cc +++ b/Helicity/Vertex/Tensor/FFTVertex.cc @@ -1,296 +1,296 @@ // -*- C++ -*- // // FFTVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FFTVertex class. // #include "FFTVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // Definition of the static class description member // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGFFTVertex("ThePEG::FFTVertex", "libThePEG.so"); void FFTVertex::Init() { static ClassDocumentation documentation ("The FFTVertex class is the implementation of" "the fermion-antifermion tensor vertices for helicity " "amplitude calculations. All such vertices should inherit" "from it."); } // function to evaluate the vertex Complex FFTVertex::evaluate(Energy2 q2,const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, const TensorWaveFunction & ten) { // set the couplings setCoupling(q2,sp.particle(),sbar.particle(),ten.particle()); // vector current LorentzPolarizationVector as = sp.wave().vectorCurrent(sbar.wave()); // momentum difference Lorentz5Momentum vdiff = sp.momentum()-sbar.momentum(); // first term LorentzPolarizationVectorE test = ten.wave().postDot(vdiff) + ten.wave().preDot(vdiff); complex term1 = as.dot(test); // trace of polarization tensor Complex trace = ten.wave().trace(); // dot products with polarization tensor // product of spinors Complex ffbar= sp.wave().scalar(sbar.wave()); // put everything together return -0.125*Complex(0.,1.)*norm()*UnitRemoval::InvE* ( term1 + 4.0*(sp.particle()->mass())*trace*ffbar ); } // member function to evaluate an off-shell spinor SpinorWaveFunction FFTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const SpinorWaveFunction & sp, const TensorWaveFunction & ten, complex mass, complex width) { // momentum of the outgoing fermion Lorentz5Momentum pout = ten.momentum()+sp.momentum(); // set the couplings setCoupling(q2,sp.particle(),out,ten.particle()); Complex ii(0.,1.); // trace of the tensor Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // mass of the fermion if(mass.real() < ZERO) mass = out->mass(); // overall factor Energy2 p2 = pout.m2(); Complex fact = 0.125*norm()*propagator(iopt,p2,out,mass,width); // compute the vector we need complex dot[4]; for(int ix=0;ix<4;++ix) { // evaluate the products we need dot[ix] =(ten(ix,3)+ten(3,ix))*(pout.e()+sp.e()); dot[ix]-=(ten(ix,0)+ten(0,ix))*(pout.x()+sp.px()); dot[ix]-=(ten(ix,1)+ten(1,ix))*(pout.y()+sp.py()); dot[ix]-=(ten(ix,2)+ten(2,ix))*(pout.z()+sp.pz()); } LorentzVector > vec(dot[0],dot[1],dot[2],dot[3]); vec -= 2.0 * trace * (pout + sp.momentum()); // combinations of the vector complex a1p2=vec.x()+ii*vec.y(); complex a1m2=vec.x()-ii*vec.y(); LorentzSpinor spt =sp .wave(); // now compute the first stage of the spinor wavefunction complex a0p3=vec.t()+vec.z(); complex a0m3=vec.t()-vec.z(); vec.setX(a0m3*spt.s3()-a1m2*spt.s4()); vec.setY(a0p3*spt.s4()-a1p2*spt.s3()); vec.setZ(a0p3*spt.s1()+a1m2*spt.s2()); vec.setT(a0m3*spt.s2()+a1p2*spt.s1()); if(mass.real()!=ZERO) { complex dot = 4.*mass*trace; vec.setX(vec.x() + dot*spt.s1()); vec.setY(vec.y() + dot*spt.s2()); vec.setZ(vec.z() + dot*spt.s3()); vec.setT(vec.t() + dot*spt.s4()); } // combinations of the momentum complex p1p2=pout.x()+ii*pout.y(); complex p1m2=pout.x()-ii*pout.y(); // finally put everything together as the spinor Complex ferm[4]; complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); - ferm[0] = UnitRemoval::InvE2 * fact*( p0m3*vec.z()-p1m2*vec.t()); - ferm[1] = UnitRemoval::InvE2 * fact*(-p1p2*vec.z()+p0p3*vec.t()); - ferm[2] = UnitRemoval::InvE2 * fact*( p0p3*vec.x()+p1m2*vec.y()); - ferm[3] = UnitRemoval::InvE2 * fact*( p1p2*vec.x()+p0m3*vec.y()); + ferm[0] = Complex(UnitRemoval::InvE2 * fact*( p0m3*vec.z()-p1m2*vec.t())); + ferm[1] = Complex(UnitRemoval::InvE2 * fact*(-p1p2*vec.z()+p0p3*vec.t())); + ferm[2] = Complex(UnitRemoval::InvE2 * fact*( p0p3*vec.x()+p1m2*vec.y())); + ferm[3] = Complex(UnitRemoval::InvE2 * fact*( p1p2*vec.x()+p0m3*vec.y())); if(mass.real()!=ZERO) { - ferm[0] += UnitRemoval::InvE2 * fact*(mass*vec.x()); - ferm[1] += UnitRemoval::InvE2 * fact*(mass*vec.y()); - ferm[2] += UnitRemoval::InvE2 * fact*(mass*vec.z()); - ferm[3] += UnitRemoval::InvE2 * fact*(mass*vec.t()); + ferm[0] += Complex(UnitRemoval::InvE2 * fact*(mass*vec.x())); + ferm[1] += Complex(UnitRemoval::InvE2 * fact*(mass*vec.y())); + ferm[2] += Complex(UnitRemoval::InvE2 * fact*(mass*vec.z())); + ferm[3] += Complex(UnitRemoval::InvE2 * fact*(mass*vec.t())); } // return the wavefunction return SpinorWaveFunction(pout,out,ferm[0],ferm[1],ferm[2],ferm[3]); } // member function to evaluate an off-shell spinor bar SpinorBarWaveFunction FFTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const SpinorBarWaveFunction & sbar, const TensorWaveFunction & ten, complex mass, complex width) { // momentum of the outgoing fermion Lorentz5Momentum pout = ten.momentum()+sbar.momentum(); // set the couplings setCoupling(q2,out,sbar.particle(),ten.particle()); Complex ii(0.,1.); // trace of the tensor Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // mass of the fermion if(mass.real() < ZERO) mass = out->mass(); // overall factor Energy2 p2 = pout.m2(); Complex fact=0.125*norm()*propagator(iopt,p2,out,mass,width); // vector complex dot[4]; for(int ix=0;ix<4;++ix) { // evaluate the products we need dot[ix] =-(ten(ix,3)+ten(3,ix))*(pout.e()+sbar.e()); dot[ix]+= (ten(ix,0)+ten(0,ix))*(pout.x()+sbar.px()); dot[ix]+= (ten(ix,1)+ten(1,ix))*(pout.y()+sbar.py()); dot[ix]+= (ten(ix,2)+ten(2,ix))*(pout.z()+sbar.pz()); } LorentzVector > vec(dot[0],dot[1],dot[2],dot[3]); vec += 2.*trace*(pout+sbar.momentum()); // combinations of the vector complex a1p2=vec.x()+ii*vec.y(); complex a1m2=vec.x()-ii*vec.y(); LorentzSpinorBar sbart=sbar.wave(); // now compute the first stage of the spinorbar wavefunction complex a0p3=vec.t()+vec.z(); complex a0m3=vec.t()-vec.z(); vec.setX(a0p3*sbart.s3()+a1p2*sbart.s4()); vec.setY(a0m3*sbart.s4()+a1m2*sbart.s3()); vec.setZ(a0m3*sbart.s1()-a1p2*sbart.s2()); vec.setT(a0p3*sbart.s2()-a1m2*sbart.s1()); if(mass.real()!=ZERO) { complex dot = 4.*mass*trace; vec.setX(vec.x() + dot*sbart.s1()); vec.setY(vec.y() + dot*sbart.s2()); vec.setZ(vec.z() + dot*sbart.s3()); vec.setT(vec.t() + dot*sbart.s4()); } // combinations of the momentum complex p1p2=pout.x()+ii*pout.y(); complex p1m2=pout.x()-ii*pout.y(); // finally put everything together as the spinor Complex ferm[4]; complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); - ferm[0] = UnitRemoval::InvE2 * fact*(-p0p3*vec.z()-p1p2*vec.t()); - ferm[1] = UnitRemoval::InvE2 * fact*(-p1m2*vec.z()-p0m3*vec.t()); - ferm[2] = UnitRemoval::InvE2 * fact*(-p0m3*vec.x()+p1p2*vec.y()); - ferm[3] = UnitRemoval::InvE2 * fact*( p1m2*vec.x()-p0p3*vec.y()); + ferm[0] = Complex(UnitRemoval::InvE2 * fact*(-p0p3*vec.z()-p1p2*vec.t())); + ferm[1] = Complex(UnitRemoval::InvE2 * fact*(-p1m2*vec.z()-p0m3*vec.t())); + ferm[2] = Complex(UnitRemoval::InvE2 * fact*(-p0m3*vec.x()+p1p2*vec.y())); + ferm[3] = Complex(UnitRemoval::InvE2 * fact*( p1m2*vec.x()-p0p3*vec.y())); if(mass.real()!=ZERO) { - ferm[0] += UnitRemoval::InvE2 * fact*mass*vec.x(); - ferm[1] += UnitRemoval::InvE2 * fact*mass*vec.y(); - ferm[2] += UnitRemoval::InvE2 * fact*mass*vec.z(); - ferm[3] += UnitRemoval::InvE2 * fact*mass*vec.t(); + ferm[0] += Complex(UnitRemoval::InvE2 * fact*mass*vec.x()); + ferm[1] += Complex(UnitRemoval::InvE2 * fact*mass*vec.y()); + ferm[2] += Complex(UnitRemoval::InvE2 * fact*mass*vec.z()); + ferm[3] += Complex(UnitRemoval::InvE2 * fact*mass*vec.t()); } // return the wavefunction return SpinorBarWaveFunction(pout,out,ferm[0],ferm[1],ferm[2],ferm[3]); } // member function to evaluate an off-shell tensor TensorWaveFunction FFTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, complex mass, complex width) { // calculating the couplings setCoupling(q2,sp.particle(),sbar.particle(),out); Complex ii(0.,1.); // momentum of the outgoing tensor Lorentz5Momentum pout = sp.momentum()+sbar.momentum(); // calculate the prefactor Energy2 p2 = pout.m2(); Complex fact = 0.0625*norm()*propagator(iopt,p2,out,mass,width); if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); // spinor vector Complex aspin[4]; LorentzSpinorBar sbart=sbar.wave(); LorentzSpinor spt =sp .wave(); aspin[3] = sbart.s1()*spt.s3()+sbart.s2()*spt.s4() +sbart.s3()*spt.s1()+sbart.s4()*spt.s2(); // spatial components are the same in both conventions aspin[0] = +sbart.s1()*spt.s4()+sbart.s2()*spt.s3() -sbart.s3()*spt.s2()-sbart.s4()*spt.s1(); aspin[1] = ii*(-sbart.s1()*spt.s4()+sbart.s2()*spt.s3() +sbart.s3()*spt.s2()-sbart.s4()*spt.s1()); aspin[2] = +sbart.s1()*spt.s3()-sbart.s2()*spt.s4() -sbart.s3()*spt.s1()+sbart.s4()*spt.s2(); // mass dependent term Complex ffbar; if(sp.particle()->mass()!=ZERO) { - ffbar = UnitRemoval::InvE * (sp.particle()->mass())* - (sp.s1()*sbar.s1()+sp.s2()*sbar.s2()+sp.s3()*sbar.s3()+sp.s4()*sbar.s4()); + ffbar = Complex(UnitRemoval::InvE * (sp.particle()->mass())* + (sp.s1()*sbar.s1()+sp.s2()*sbar.s2()+sp.s3()*sbar.s3()+sp.s4()*sbar.s4())); } else { ffbar = 0.; } // dot products for the calculation complex dotka = +aspin[3]*pout.e()-aspin[0]*pout.x() -aspin[1]*pout.y()-aspin[2]*pout.z(); complex dot12a = +aspin[3]*(sp.e() -sbar.e() )-aspin[0]*(sp.px()-sbar.px()) -aspin[1]*(sp.py()-sbar.py())-aspin[2]*(sp.pz()-sbar.pz()); complex diff=sp.m2()-sbar.m2(); complex dotkam=dotka/mass2; Complex diffm =diff/mass2; Complex p2m = p2/mass2; // construct the vectors for the first two terms Complex veca[4],vecb[4]; veca[0] = aspin[0]-UnitRemoval::InvE2*dotka*pout.x(); - vecb[0] = UnitRemoval::InvE*(sp.px()-sbar.px()-diffm*pout.x()); + vecb[0] = Complex(UnitRemoval::InvE*(sp.px()-sbar.px()-diffm*pout.x())); veca[1] = aspin[1]-UnitRemoval::InvE2*dotka*pout.y(); - vecb[1] = UnitRemoval::InvE*(sp.py()-sbar.py()-diffm*pout.y()); + vecb[1] = Complex(UnitRemoval::InvE*(sp.py()-sbar.py()-diffm*pout.y())); veca[2] = aspin[2]-UnitRemoval::InvE2*dotka*pout.z(); - vecb[2] = UnitRemoval::InvE*(sp.pz()-sbar.pz()-diffm*pout.z()); + vecb[2] = Complex(UnitRemoval::InvE*(sp.pz()-sbar.pz()-diffm*pout.z())); veca[3] = aspin[3]-UnitRemoval::InvE2*dotka*pout.e(); - vecb[3] = UnitRemoval::InvE*(sp.e()-sbar.e()-diffm*pout.e()); + vecb[3] = Complex(UnitRemoval::InvE*(sp.e()-sbar.e()-diffm*pout.e())); // coefficients fr hte second two terms Complex temp = UnitRemoval::InvE*(p2m*dot12a-dotkam*diff); Complex coeff1 = -4./3.*(2.*ffbar*(1.-p2m) + temp); - temp = UnitRemoval::InvE*(-3.*dot12a+2.*p2m*dot12a+diffm*dotka); + temp = Complex(UnitRemoval::InvE*(-3.*dot12a+2.*p2m*dot12a+diffm*dotka)); Complex coeff2 = -4./3./mass2*( 4.*ffbar*(1.-p2m) + temp)*UnitRemoval::E2; // construct the tensor Complex ten[4][4]; const complex pout_tmp[4] = {pout.x(), pout.y(), pout.z(), pout.e()}; for(int ix=0;ix<4;++ix) { for(int iy=0;iy<4;++iy) { Complex temp = coeff2*pout_tmp[ix]*pout_tmp[iy]*UnitRemoval::InvE2; ten[ix][iy] = 2.*(veca[ix]*vecb[iy]+veca[iy]*vecb[ix]) + temp; } } ten[0][0]=ten[0][0]-coeff1; ten[1][1]=ten[1][1]-coeff1; ten[2][2]=ten[2][2]-coeff1; ten[3][3]=ten[3][3]+coeff1; // multiply by final prefactor for(int ix=0;ix<4;++ix) { for(int iy=0;iy<4;++iy) { ten[ix][iy] = fact*ten[ix][iy]; } } // return the wavefunction return TensorWaveFunction(pout,out, ten[0][0],ten[0][1],ten[0][2],ten[0][3], ten[1][0],ten[1][1],ten[1][2],ten[1][3], ten[2][0],ten[2][1],ten[2][2],ten[2][3], ten[3][0],ten[3][1],ten[3][2],ten[3][3]); } diff --git a/Helicity/Vertex/Tensor/SSTVertex.cc b/Helicity/Vertex/Tensor/SSTVertex.cc --- a/Helicity/Vertex/Tensor/SSTVertex.cc +++ b/Helicity/Vertex/Tensor/SSTVertex.cc @@ -1,161 +1,161 @@ // -*- C++ -*- // // SSTVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the SSTVertex class. // #include "SSTVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGSSTVertex("ThePEG::SSTVertex", "libThePEG.so"); void SSTVertex::Init() { static ClassDocumentation documentation ("The SSTVertex class is the implementation of the " "helicity amplitude calculation for the scalar-scalar-tensor vertex" ", all such vertices should inherit from it"); } // evaluate the vertex Complex SSTVertex::evaluate(Energy2 q2, const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2, const TensorWaveFunction & ten) { // obtain the coupling setCoupling(q2,sca1.particle(),sca2.particle(),ten.particle()); Complex ii(0.,1.); // evaluate the trace of the tensor Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // dot product of the two momenta Energy2 dot = sca1.momentum()*sca2.momentum(); Energy mass = sca1.particle()->mass(); // second term complex second = +2.*ten.tt()*sca1.e()*sca2.e() +2.*ten.xx()*sca1.px()*sca2.px() +2.*ten.yy()*sca1.py()*sca2.py()+2.*ten.zz()*sca1.pz()*sca2.pz() -(ten.tx()+ten.xt())*(sca1.e()*sca2.px() +sca1.px()*sca2.e()) -(ten.ty()+ten.yt())*(sca1.e()*sca2.py() +sca1.py()*sca2.e()) -(ten.tz()+ten.zt())*(sca1.e()*sca2.pz() +sca1.pz()*sca2.e()) +(ten.xy()+ten.yx())*(sca1.py()*sca2.px()+sca1.px()*sca2.py()) +(ten.xz()+ten.zx())*(sca1.pz()*sca2.px()+sca1.px()*sca2.pz()) +(ten.yz()+ten.zy())*(sca1.pz()*sca2.py()+sca1.py()*sca2.pz()); // return the answer return -0.5*ii*norm()*UnitRemoval::InvE2* (trace*(mass*mass-dot)+second)*sca1.wave()*sca2.wave(); } // off-shell tensor TensorWaveFunction SSTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const ScalarWaveFunction & sca1, const ScalarWaveFunction & sca2, complex mass, complex width) { // obtain the coupling setCoupling(q2,sca1.particle(),sca2.particle(),out); // array for the tensor Complex ten[4][4]; // calculate the outgoing momentum Lorentz5Momentum pout = sca1.momentum()+sca2.momentum(); // prefactor if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); Energy2 p2 = pout.m2(); Complex fact = 0.25*norm()*sca1.wave()*sca2.wave()* propagator(iopt,p2,out,mass,width); // dot products we need Energy2 dot12 = sca1.momentum()*sca2.momentum(); Energy2 dot1 = sca1.momentum()*pout; Energy2 dot2 = pout*sca2.momentum(); // the vectors that we need for the tensor LorentzPolarizationVectorE vec1,vec2; Complex a,b; Energy2 mphi2 = sqr(sca1.particle()->mass()); // massive case if(mass.real()!=ZERO) { Complex norm1=dot1/mass2; Complex norm2=dot2/mass2; a = Complex(UnitRemoval::InvE2 * ((mphi2+dot12)*(Complex(2.*p2/mass2)-5.) +4.*(dot12-dot1*dot2/mass2)))/3.; b = -(-(mphi2+dot12)*(2.+p2/mass2)+4.*(dot12-dot1*(dot2/mass2)))/3./mass2; vec1 = LorentzPolarizationVectorE(sca1.momentum()) - LorentzPolarizationVectorE(norm1 * pout); vec2 = LorentzPolarizationVectorE(sca2.momentum()) - LorentzPolarizationVectorE(norm2 * pout); } // massless case else { a = UnitRemoval::InvE2 * (-5.*(mphi2+dot12)+4.*dot12)/3.; b = 0.; vec1 = sca1.momentum(); vec2 = sca2.momentum(); } // calculate the wavefunction complex vec1_tmp[4] = {vec1.x(), vec1.y(), vec1.z(), vec1.t()}; complex vec2_tmp[4] = {vec2.x(), vec2.y(), vec2.z(), vec2.t()}; complex pout_tmp[4] = {pout.x(), pout.y(), pout.z(), pout.t()}; for(int ix=0;ix<4;++ix) { for(int iy=0;iy<4;++iy) { complex temp = -2.*( vec1_tmp[ix]*vec2_tmp[iy] + vec1_tmp[ix]*vec2_tmp[iy]) -b*pout_tmp[ix]*pout_tmp[iy]; - ten[ix][iy]= UnitRemoval::InvE2 * temp; + ten[ix][iy]= Complex(UnitRemoval::InvE2 * temp); } } ten[3][3]=ten[3][3]-a; for(int ix=0;ix<3;++ix){ten[ix][ix]=ten[ix][ix]+a;} // prefactor for(int ix=0;ix<4;++ix){for(int iy=0;iy<4;++iy){ten[ix][iy]=fact*ten[ix][iy];}} // return the wavefunction return TensorWaveFunction(pout,out, ten[0][0],ten[0][1],ten[0][2],ten[0][3], ten[1][0],ten[1][1],ten[1][2],ten[1][3], ten[2][0],ten[2][1],ten[2][2],ten[2][3], ten[3][0],ten[3][1],ten[3][2],ten[3][3]); } // off-shell scalar ScalarWaveFunction SSTVertex::evaluate(Energy2 q2,int iopt, tcPDPtr out, const ScalarWaveFunction & sca, const TensorWaveFunction & ten, complex mass, complex width) { // obtain the coupling setCoupling(q2,sca.particle(),out,ten.particle()); // calculate the outgoing momentum Lorentz5Momentum pout = sca.momentum()+ten.momentum(); // prefactors if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); Energy2 p2 = pout.m2(); Complex fact = 0.5*norm()*sca.wave()*propagator(iopt,p2,out,mass,width); // trace of the tensor Complex trace1 = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // dot product of the two momenta Energy2 dot = sca.momentum()*pout; // first term complex trace = trace1*(mass2-dot); // second term complex second = +2.*ten.tt()*sca.e()*pout.e() +2.*ten.xx()*sca.px()*pout.x() +2.*ten.yy()*sca.py()*pout.y()+2.*ten.zz()*sca.pz()*pout.z() -(ten.tx()+ten.xt())*( sca.e()*pout.x()+sca.px()*pout.e()) -(ten.ty()+ten.yt())*( sca.e()*pout.y()+sca.py()*pout.e()) -(ten.tz()+ten.zt())*( sca.e()*pout.z()+sca.pz()*pout.e()) +(ten.xy()+ten.yx())*(sca.py()*pout.x()+sca.px()*pout.y()) +(ten.xz()+ten.zx())*(sca.pz()*pout.x()+sca.px()*pout.z()) +(ten.yz()+ten.zy())*(sca.py()*pout.z()+sca.pz()*pout.y()); // put it all together second = fact*(trace+second); Complex result = second * UnitRemoval::InvE2; return ScalarWaveFunction(pout,out,result); } diff --git a/Helicity/Vertex/Tensor/VVTVertex.cc b/Helicity/Vertex/Tensor/VVTVertex.cc --- a/Helicity/Vertex/Tensor/VVTVertex.cc +++ b/Helicity/Vertex/Tensor/VVTVertex.cc @@ -1,277 +1,277 @@ // -*- C++ -*- // // VVTVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the VVTVertex class. // #include "VVTVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGVVTVertex("ThePEG::VVTVertex", "libThePEG.so"); void VVTVertex::Init() { static ClassDocumentation documentation ("The VVTVertex class is the implementation of" "the vector-vector tensor vertices for helicity " "amplitude calculations. All such vertices should inherit" "from it."); } // function to evaluate the vertex Complex VVTVertex::evaluate(Energy2 q2, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const TensorWaveFunction & ten, Energy vmass) { // set the couplings setCoupling(q2,vec1.particle(),vec2.particle(),ten.particle()); // mass of the vector if(vmassmass(); // mass+k1.k2 Energy2 mdot = vec1.momentum()*vec2.momentum(); if(vmass!=ZERO) mdot += sqr(vmass); // dot product of wavefunctions and momenta Complex dotv1v2 = vec1.wave().dot(vec2.wave()); complex dotk1v2 = vec1.momentum()*vec2.wave(); complex dotk2v1 = vec1.wave()*vec2.momentum(); // dot product of wavefunctions and momenta with the tensor LorentzPolarizationVector tv1 = ten.wave().postDot(vec1.wave()); LorentzPolarizationVector tv2 = ten.wave().postDot(vec2.wave()); LorentzPolarizationVectorE tk1 = ten.wave().postDot(vec1.momentum()); LorentzPolarizationVectorE tk2 = ten.wave().postDot(vec2.momentum()); Complex tenv1v2 = vec1.wave ().dot(tv2) + vec2.wave ().dot(tv1); complex tenk1v2 = vec1.momentum().dot(tv2) + vec2.wave ().dot(tk1); complex tenk2v1 = vec2.momentum().dot(tv1) + vec1.wave ().dot(tk2); complex tenk1k2 = vec2.momentum().dot(tk1) + vec1.momentum().dot(tk2); // trace of the tensor Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // evaluate the vertex return -0.5*Complex(0.,1.)*norm()*UnitRemoval::InvE2 * (trace*(dotk1v2*dotk2v1-dotv1v2*mdot) +mdot*tenv1v2-dotk2v1*tenk1v2 -dotk1v2*tenk2v1+dotv1v2*tenk1k2); } // evaluate an off-shell vector VectorWaveFunction VVTVertex::evaluate(Energy2 q2, int iopt, tcPDPtr out, const VectorWaveFunction & vec, const TensorWaveFunction & ten, complex mass, complex width) { // evaluate the couplings setCoupling(q2,vec.particle(),out,ten.particle()); // outgoing momentum Lorentz5Momentum pout = ten.momentum()+vec.momentum(); // normalisation factor if(mass.real() < ZERO) mass = out->mass(); complex mass2 = sqr(mass); Energy2 p2 = pout.m2(); Complex fact = -0.5*norm()*propagator(iopt,p2,out,mass,width); // dot product of wavefunctions and momenta complex dotk1k2 = vec.momentum()*pout; complex dotk2v1 = vec.wave() *pout; // mass-k1.k2 complex mdot = -dotk1k2; if(mass.real()!=ZERO) mdot += mass2; // components of the tensor Complex tentx = ten.tx()+ten.xt(); Complex tenty = ten.ty()+ten.yt(); Complex tentz = ten.tz()+ten.zt(); Complex tenxy = ten.xy()+ten.yx(); Complex tenxz = ten.xz()+ten.zx(); Complex tenyz = ten.yz()+ten.zy(); // dot product of momenta and polarization vectors with the tensor complex tenk2v1 = 2.*(+ten.tt()*vec.t()*pout.e() +ten.xx()*vec.x()*pout.x() +ten.yy()*vec.y()*pout.y()+ten.zz()*vec.z()*pout.z()) -tentx*(vec.t()*pout.x()+vec.x()*pout.e()) -tenty*(vec.t()*pout.y()+vec.y()*pout.e()) -tentz*(vec.t()*pout.z()+vec.z()*pout.e()) +tenxy*(vec.x()*pout.y()+vec.y()*pout.x()) +tenxz*(vec.x()*pout.z()+vec.z()*pout.x()) +tenyz*(vec.y()*pout.z()+vec.z()*pout.y()); complex tenk1k2 = 2.*(+ten.tt()*vec.e()*pout.e() +ten.xx()*vec.px()*pout.x() +ten.yy()*vec.py()*pout.y()+ten.zz()*vec.pz()*pout.z()) -tentx*(vec.e()*pout.x() +vec.px()*pout.e()) -tenty*(vec.e()*pout.y() +vec.py()*pout.e()) -tentz*(vec.e()*pout.z() +vec.pz()*pout.e()) +tenxy*(vec.px()*pout.y()+vec.py()*pout.x()) +tenxz*(vec.px()*pout.z()+vec.pz()*pout.x()) +tenyz*(vec.py()*pout.z()+vec.pz()*pout.y()); // trace of the tensor Complex trace = ten.tt()-ten.xx()-ten.yy()-ten.zz(); // compute the vector Complex vec1[4]; - vec1[0] = UnitRemoval::InvE2* + vec1[0] = Complex(UnitRemoval::InvE2* (mdot*(+ tentx* vec.t()-2.*ten.xx()*vec.x() - tenxy* vec.y()- tenxz* vec.z()-trace*vec.x()) +(tenk2v1-trace*dotk2v1)*vec.px()-tenk1k2*vec.x() +dotk2v1*(+ tentx* vec.e() -2.*ten.xx()*vec.px() - - tenxy* vec.py()- tenxz* vec.pz())); - vec1[1] = UnitRemoval::InvE2*( + - tenxy* vec.py()- tenxz* vec.pz()))); + vec1[1] = Complex(UnitRemoval::InvE2*( mdot * (+ tenty* vec.t() - tenxy* vec.x() -2.*ten.yy()*vec.y() - tenyz* vec.z()-trace*vec.y() ) +(tenk2v1-trace*dotk2v1)*vec.py() -tenk1k2*vec.y() +dotk2v1*(+ tenty* vec.e() - tenxy* vec.px() -2.*ten.yy()*vec.py() - - tenyz* vec.pz())); + - tenyz* vec.pz()))); - vec1[2] = UnitRemoval::InvE2* + vec1[2] = Complex(UnitRemoval::InvE2* (mdot* (+ tentz* vec.t()- tenxz* vec.x() - tenyz* vec.y()-2.*ten.zz()*vec.z()-trace*vec.z()) +(tenk2v1-trace*dotk2v1)*vec.pz()-tenk1k2*vec.z() +dotk2v1*(+ tentz* vec.e() - tenxz *vec.px() - - tenyz* vec.py()-2.*ten.zz()*vec.pz())); - vec1[3] = UnitRemoval::InvE2* + - tenyz* vec.py()-2.*ten.zz()*vec.pz()))); + vec1[3] = Complex(UnitRemoval::InvE2* (mdot*(+2.*ten.tt()*vec.t()- tentx* vec.x() - tenty* vec.y()- tentz* vec.z()-trace*vec.t()) +(tenk2v1-trace*dotk2v1)*vec.e() -tenk1k2*vec.t() +dotk2v1*(+2.*ten.tt()*vec.e() - tentx* vec.px() - - tenty* vec.py()- tentz* vec.pz())); + - tenty* vec.py()- tentz* vec.pz()))); // now add the piece for massive bosons if(mass.real()!=ZERO) { // DGRELL unit problem? Complex dot = tenk2v1 * UnitRemoval::InvE - dotk1k2 * trace * UnitRemoval::InvE2; - vec1[0] -= dot*pout.x() * UnitRemoval::InvE; - vec1[1] -= dot*pout.y() * UnitRemoval::InvE; - vec1[2] -= dot*pout.z() * UnitRemoval::InvE; - vec1[3] -= dot*pout.e() * UnitRemoval::InvE; + vec1[0] -= Complex(dot*pout.x() * UnitRemoval::InvE); + vec1[1] -= Complex(dot*pout.y() * UnitRemoval::InvE); + vec1[2] -= Complex(dot*pout.z() * UnitRemoval::InvE); + vec1[3] -= Complex(dot*pout.e() * UnitRemoval::InvE); } // return the VectorWaveFunction for(int ix=0;ix<4;++ix){vec1[ix]=vec1[ix]*fact;} return VectorWaveFunction(pout,out,vec1[0],vec1[1],vec1[2],vec1[3]); } // off-shell tensor TensorWaveFunction VVTVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, Energy vmass, complex tmass, complex width) { // coupling setCoupling(q2,vec1.particle(),vec2.particle(),out); // momenta of the outgoing tensor // outgoing momentum Lorentz5Momentum pout= vec1.momentum()+vec2.momentum(); // overall normalisation if(tmass.real() < ZERO) tmass = out->mass(); complex tmass2 = sqr(tmass); if(vmassmass(); Energy2 vmass2 = sqr(vmass); Energy2 p2 = pout.m2(); Complex fact = 0.25*norm()*propagator(iopt,p2,out,tmass,width); // dot products we need to construct the tensor complex dotk1k2 = vec1.momentum()*vec2.momentum(); complex dotv1k2 = vec1.wave()*vec2.momentum(); Complex dotv1v2 = vec1.wave().dot(vec2.wave()); complex dotk1v2 = vec1.momentum()*vec2.wave(); complex dotkk1 = vec1.momentum()*pout; complex dotkk2 = vec2.momentum()*pout; complex dotkv1 = vec1.wave()*pout; complex dotkv2 = vec2.wave()*pout; // dot product ma^2+k1.k2 complex mdot = vmass2+dotk1k2; // vectors to help construct the tensor Complex vecv1[4],vecv2[4]; complex veck1[4],veck2[4]; complex tmass2inv(ZERO); if(tmass.real()> ZERO) tmass2inv = 1./tmass2; vecv1[0]=vec1.x() -pout.x()*dotkv1*tmass2inv; vecv2[0]=vec2.x() -pout.x()*dotkv2*tmass2inv; veck1[0]=vec1.px()-pout.x()*dotkk1*tmass2inv; veck2[0]=vec2.px()-pout.x()*dotkk2*tmass2inv; vecv1[1]=vec1.y() -pout.y()*dotkv1*tmass2inv; vecv2[1]=vec2.y() -pout.y()*dotkv2*tmass2inv; veck1[1]=vec1.py()-pout.y()*dotkk1*tmass2inv; veck2[1]=vec2.py()-pout.y()*dotkk2*tmass2inv; vecv1[2]=vec1.z() -pout.z()*dotkv1*tmass2inv; vecv2[2]=vec2.z() -pout.z()*dotkv2*tmass2inv; veck1[2]=vec1.pz()-pout.z()*dotkk1*tmass2inv; veck2[2]=vec2.pz()-pout.z()*dotkk2*tmass2inv; vecv1[3]=vec1.t() -pout.e()*dotkv1*tmass2inv; vecv2[3]=vec2.t() -pout.e()*dotkv2*tmass2inv; veck1[3]=vec1.e() -pout.e()*dotkk1*tmass2inv; veck2[3]=vec2.e() -pout.e()*dotkk2*tmass2inv; // coefficient of g(nu,mu)-k^muk^nu/m^2 Complex omp2 = 1.-Complex(p2*tmass2inv); Complex coeff1 = UnitRemoval::InvE2 * ( +4./3.*mdot*(-2.*dotv1v2+Complex((dotkv1*dotkv2+p2*dotv1v2)*tmass2inv)) +4./3.*(dotv1k2*(dotk1v2-dotkk1*dotkv2*tmass2inv) +dotk1v2*(dotv1k2-dotkk2*dotkv1*tmass2inv) -dotv1v2*(dotk1k2-dotkk1*dotkk2*tmass2inv) +omp2*dotk1v2*dotv1k2) ); // coefficient of g(nu,mu) Complex coeff2 = UnitRemoval::InvE2 * ( 2.*mdot*omp2*dotv1v2 -2.*omp2*dotk1v2*dotv1k2 ); // construct the tensor Complex ten[4][4]; const Energy pout_tmp[4] = {pout.x(), pout.y(), pout.z(), pout.e()}; for(int ix=0;ix<4;++ix) { for(int iy=0;iy<4;++iy) { complex temp; temp = 2.*mdot* (vecv1[ix]*vecv2[iy]+vecv1[iy]*vecv2[ix]); temp -= 2.*dotv1k2*(veck1[ix]*vecv2[iy]+veck1[iy]*vecv2[ix]); temp -= 2.*dotk1v2*(veck2[ix]*vecv1[iy]+veck2[iy]*vecv1[ix]); temp += 2.*dotv1v2*(veck1[ix]*veck2[iy]+veck1[iy]*veck2[ix]); - ten[ix][iy] = UnitRemoval::InvE2 * temp + ten[ix][iy] = Complex(UnitRemoval::InvE2 * temp) -coeff1*tmass2inv*pout_tmp[ix]*pout_tmp[iy]; } } // add the g(mu,nu) term ten[0][0] = ten[0][0]-(coeff1+coeff2); ten[1][1] = ten[1][1]-(coeff1+coeff2); ten[2][2] = ten[2][2]-(coeff1+coeff2); ten[3][3] = ten[3][3]+(coeff1+coeff2); // overall coefficent for(int ix=0;ix<4;++ix) { for(int iy=0;iy<4;++iy) { ten[ix][iy] = fact*ten[ix][iy]; } } // return the wavefunction return TensorWaveFunction(pout,out, ten[0][0],ten[0][1],ten[0][2],ten[0][3], ten[1][0],ten[1][1],ten[1][2],ten[1][3], ten[2][0],ten[2][1],ten[2][2],ten[2][3], ten[3][0],ten[3][1],ten[3][2],ten[3][3]); } diff --git a/Helicity/Vertex/Vector/FFVVertex.cc b/Helicity/Vertex/Vector/FFVVertex.cc --- a/Helicity/Vertex/Vector/FFVVertex.cc +++ b/Helicity/Vertex/Vector/FFVVertex.cc @@ -1,475 +1,475 @@ // -*- C++ -*- // // FFVVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the FFVVertex class. // #include "FFVVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // Definition of the static class description member // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGFFVVertex("ThePEG::FFVVertex", "libThePEG.so"); void FFVVertex::Init() { static ClassDocumentation documentation ("The FFVVertex class implements the helicity amplitude" "calculations for a fermion-fantifermion gauge boson vertex. Any " "implementation of such a vertex should inherit from in and implement" " the virtual setCoupling member to calculate the coupling"); } // evalulate the full vertex Complex FFVVertex::evaluate(Energy2 q2, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, const VectorWaveFunction & vec) { // first calculate the couplings if(kinematics()) calculateKinematics(sp.momentum(),sbar.momentum(),vec.momentum()); setCoupling(q2,sp.particle(),sbar.particle(),vec.particle()); Complex ii(0.,1.); // useful combinations of the polarization vector components Complex e0p3=vec.t()+vec.z(); Complex e0m3=vec.t()-vec.z(); Complex e1p2=vec.x()+ii*vec.y(); Complex e1m2=vec.x()-ii*vec.y(); Complex vertex(0.); if(_left!=0.) { vertex += _left*(+sbar.s3()*(sp.s1()*e0p3+sp.s2()*e1m2) +sbar.s4()*(sp.s1()*e1p2+sp.s2()*e0m3)); } // then the right piece (often not needed eg W vertex) if(_right!=0.) { vertex += _right*(+sbar.s1()*(sp.s3()*e0m3-sp.s4()*e1m2) -sbar.s2()*(sp.s3()*e1p2-sp.s4()*e0p3)); } vertex*=ii; return vertex*norm(); } // evaluate an off-shell spinor SpinorWaveFunction FFVVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out, const SpinorWaveFunction & sp, const VectorWaveFunction &vec, complex mass, complex width) { // extract the pointers to the particle data objects tcPDPtr Psp=sp.particle(); tcPDPtr Pvec=vec.particle(); // work out the momentum of the off-shell particle Lorentz5Momentum pout = sp.momentum()+vec.momentum(); // first calculate the couplings if(kinematics()) calculateKinematics(sp.momentum(),pout,vec.momentum()); setCoupling(q2,Psp,out,Pvec); Complex ii(0.,1.); // now evaluate the contribution // polarization components Complex e0p3 = vec.t() + vec.z(); Complex e0m3 = vec.t() - vec.z(); Complex e1p2 = vec.x()+ii*vec.y(); Complex e1m2 = vec.x()-ii*vec.y(); // overall factor Energy2 p2 = pout.m2(); Complex fact=-normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = iopt==5 ? ZERO : out->mass(); complex p1p2 = pout.x()+ii*pout.y(); complex p1m2 = pout.x()-ii*pout.y(); // complex nos for for the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); LorentzSpinor spt =sp .wave(); complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); // left piece if(_left!=0.) { Complex a3=_left*fact*( spt.s1()*e0p3+spt.s2()*e1m2); Complex a4=_left*fact*( spt.s1()*e1p2+spt.s2()*e0m3); - s1 +=UnitRemoval::InvE * (p0m3*a3-p1m2*a4); - s2 +=UnitRemoval::InvE * (-p1p2*a3+p0p3*a4); - s3 +=UnitRemoval::InvE * a3*mass; - s4 +=UnitRemoval::InvE * a4*mass; + s1 += Complex(UnitRemoval::InvE * (p0m3*a3-p1m2*a4)); + s2 += Complex(UnitRemoval::InvE * (-p1p2*a3+p0p3*a4)); + s3 += Complex(UnitRemoval::InvE * a3*mass); + s4 += Complex(UnitRemoval::InvE * a4*mass); } // right piece if(_right!=0.) { Complex a1=_right*fact*( spt.s3()*e0m3-spt.s4()*e1m2); Complex a2=_right*fact*(-spt.s3()*e1p2+spt.s4()*e0p3); - s1 +=UnitRemoval::InvE * a1*mass; - s2 +=UnitRemoval::InvE * a2*mass; - s3 +=UnitRemoval::InvE * (p0p3*a1+p1m2*a2); - s4 +=UnitRemoval::InvE * (p1p2*a1+p0m3*a2); + s1 += Complex(UnitRemoval::InvE * a1*mass); + s2 += Complex(UnitRemoval::InvE * a2*mass); + s3 += Complex(UnitRemoval::InvE * (p0p3*a1+p1m2*a2)); + s4 += Complex(UnitRemoval::InvE * (p1p2*a1+p0m3*a2)); } // return the wavefunction return SpinorWaveFunction(pout,out,s1,s2,s3,s4); } // evaluate an off-shell SpinorBar SpinorBarWaveFunction FFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const SpinorBarWaveFunction & sbar, const VectorWaveFunction & vec, complex mass, complex width) { // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+vec.momentum(); // first calculate the couplings if(kinematics()) calculateKinematics(pout,sbar.momentum(),vec.momentum()); setCoupling(q2,out,sbar.particle(),vec.particle()); Complex ii(0.,1.); // now evaluate the contribution // polarization components Complex e0p3=vec.t() + vec.z(); Complex e0m3=vec.t() - vec.z(); Complex e1p2=vec.x()+ii*vec.y(); Complex e1m2=vec.x()-ii*vec.y(); // overall factor Energy2 p2 = pout.m2(); Complex fact=-normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); complex p1p2=pout.x()+ii*pout.y(); complex p1m2=pout.x()-ii*pout.y(); // complex numbers for the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); LorentzSpinorBar sbart=sbar.wave(); Energy p0p3=pout.e() + pout.z(); Energy p0m3=pout.e() - pout.z(); // left piece if(_left!=0.) { Complex a1 = _left*fact*( sbart.s3()*e0p3+sbart.s4()*e1p2); Complex a2 = _left*fact*( sbart.s3()*e1m2+sbart.s4()*e0m3); - s1 += UnitRemoval::InvE*a1*mass; - s2 += UnitRemoval::InvE*a2*mass; - s3 += UnitRemoval::InvE*(-p0m3*a1+p1p2*a2); - s4 += UnitRemoval::InvE*(+p1m2*a1-p0p3*a2); + s1 += Complex(UnitRemoval::InvE*a1*mass); + s2 += Complex(UnitRemoval::InvE*a2*mass); + s3 += Complex(UnitRemoval::InvE*(-p0m3*a1+p1p2*a2)); + s4 += Complex(UnitRemoval::InvE*(+p1m2*a1-p0p3*a2)); } // right piece if(_right!=0.) { Complex a3 = _right*fact*( sbart.s1()*e0m3-sbart.s2()*e1p2); Complex a4 = _right*fact*(-sbart.s1()*e1m2+sbart.s2()*e0p3); - s1 += UnitRemoval::InvE*(-p0p3*a3-p1p2*a4); - s2 += UnitRemoval::InvE*(-p1m2*a3-p0m3*a4); - s3 += UnitRemoval::InvE*a3*mass; - s4 += UnitRemoval::InvE*a4*mass; + s1 += Complex(UnitRemoval::InvE*(-p0p3*a3-p1p2*a4)); + s2 += Complex(UnitRemoval::InvE*(-p1m2*a3-p0m3*a4)); + s3 += Complex(UnitRemoval::InvE*a3*mass); + s4 += Complex(UnitRemoval::InvE*a4*mass); } return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4); } // off-shell vector VectorWaveFunction FFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, complex mass, complex width) { // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+sp.momentum(); // first calculate the couplings if(kinematics()) calculateKinematics(sp.momentum(),sbar.momentum(),pout); setCoupling(q2,sp.particle(),sbar.particle(),out); Complex ii(0.,1.); // overall factor Energy2 p2 = pout.m2(); Complex fact = normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); complex mass2 = sqr(mass); // the vector for the fermion-antifermion Complex vec[4]; // left coupling if(_left!=0.) { vec[0] = -_left*(sbar.s3()*sp.s2()+sbar.s4()*sp.s1()); vec[1] = ii*_left*(sbar.s3()*sp.s2()-sbar.s4()*sp.s1()); vec[2] = -_left*(sbar.s3()*sp.s1()-sbar.s4()*sp.s2()); vec[3] = _left*(sbar.s3()*sp.s1()+sbar.s4()*sp.s2()); } // right coupling if(_right!=0.) { vec[0] += +_right*(sbar.s1()*sp.s4()+sbar.s2()*sp.s3()); vec[1] += -ii*_right*(sbar.s1()*sp.s4()-sbar.s2()*sp.s3()); vec[2] += +_right*(sbar.s1()*sp.s3()-sbar.s2()*sp.s4()); vec[3] += +_right*(sbar.s1()*sp.s3()+sbar.s2()*sp.s4()); } // massless boson if(mass.real()==ZERO) { for(int ix=0;ix<4;++ix){vec[ix]*=fact;} } // massive boson else { complex dot = ( pout.e() *vec[3] -pout.x()*vec[0] -pout.y()*vec[1] -pout.z()*vec[2])/mass2; vec[0]=fact*(vec[0]-dot*pout.x()); vec[1]=fact*(vec[1]-dot*pout.y()); vec[2]=fact*(vec[2]-dot*pout.z()); vec[3]=fact*(vec[3]-dot*pout.e()); } return VectorWaveFunction(pout,out,vec[0],vec[1],vec[2],vec[3]); } SpinorWaveFunction FFVVertex::evaluateSmall(Energy2 q2,int iopt, tcPDPtr out, const SpinorWaveFunction & sp, const VectorWaveFunction & vec, unsigned int fhel, unsigned int vhel, double ctheta, double phi, double stheta, bool includeEikonal, SmallAngleDirection direction, Energy mass, Energy) { assert(fhel <= 1); assert( vhel == 0 || vhel == 2 ); SpinorWaveFunction output; // first calculate the couplings setCoupling(q2,sp.particle(),out,vec.particle()); Complex ii(0.,1.); if(mass < ZERO) mass = iopt==5 ? ZERO : out->mass(); Lorentz5Momentum pout = sp.momentum()+vec.momentum(); assert(sp.direction()!=intermediate); // helicity of the boson double lam = double(vhel)-1.; // energy of the boson Energy Eg = abs(vec.momentum().e()); // energy of the fermion Energy Ef = abs(sp .momentum().e()); // energy fraction of photon double x = Eg/Ef; // velocity of the fermon double beta = sqrt(1.-sqr(mass/Ef)); // dimensionless versions of the variables double dm = mass*UnitRemoval::InvE; double dEf = Ef*UnitRemoval::InvE; double rE = sqrt(.5*dEf); // calculation of propagator accurate as beta->1 and theta -> 0 Energy2 dot = 2.*Ef*Eg*(sqr(mass/Ef)/(1.+beta)*ctheta + sqr(stheta)/(1.+ctheta) ); Complex fact= norm()*(0.5*left()+0.5*right())*UnitRemoval::E2/dot; // phase factor Complex ephig = cos(phi )+ii*sin(phi ); // calculation of the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); // u-type spinor if(sp.wave().Type()==SpinorType::u) { // fermion along +z if(direction==PostiveZDirection) { if(fhel==0) { s1 = +x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s2 =-rE*dEf*sqrt(1.+beta)*stheta* (+x*(1.+lam)-(includeEikonal ? 2.*beta*lam : 0. )); s3 = +x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta); s4 = +rE*dm/sqrt(1.+beta)*stheta* (-x*(1.-lam)+(includeEikonal ? 2.*beta*lam : 0. )); } else if(fhel==1) { s1 = +rE*dm/sqrt(1.+beta)*stheta* (+x*(1.+lam)+(includeEikonal ? 2.*beta*lam : 0. )); s2 = -x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta); s3 = -rE*dEf*sqrt(1.+beta)*stheta* (-x*(1.-lam)-(includeEikonal ? 2.*beta*lam : 0. )); s4 = -x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta); } } // fermion along -z else { if(fhel==0) { s1 = -rE*dEf*sqrt(1.+beta)*stheta* (+x*(1.+lam)-(includeEikonal ? 2.*beta*lam : 0. )); s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s3 = rE*dm/sqrt(1.+beta)*stheta* (-x*(1.-lam)+(includeEikonal ? 2.*beta*lam : 0. )); s4 = -x*rE*dm/ephig*(-1.+lam)*(1.+ctheta)/sqrt(1.+beta); } else if(fhel==1) { s1 =-x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta); s2 =-rE*dm/sqrt(1.+beta)*stheta* ( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. )); s3 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta); s4 =-rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); } } } // v-type spinor else if(sp.wave().Type()==SpinorType::v) { // fermion along +z if(direction==PostiveZDirection) { if(fhel==0) { s1 = rE*dm/sqrt(1.+beta)*stheta* (-x*(1.+lam) + ( includeEikonal ? 2.*beta*lam : 0. )); s2 = x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta); s3 = rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) - ( includeEikonal ? 2.*beta*lam : 0. )); s4 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta); } else if(fhel==1) { s1 = x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s2 =-rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0.)); s3 = x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta); s4 = rE*dm/sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); } } // fermion aling -z else { if(fhel==0) { s1 = x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta); s2 = rE*dm/sqrt(1.+beta)*stheta* ( x*(1.+lam) - ( includeEikonal ? 2.*beta*lam : 0. )); s3 = x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta) / (1.+ctheta); s4 =-rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) - ( includeEikonal ? 2.*beta*lam : 0. )); } else if(fhel==1) { s1 =-rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.+lam) + ( includeEikonal ? 2.*beta*lam : 0. )); s2 =-x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s3 = rE*dm/sqrt(1.+beta)*stheta* ( x*(1.-lam) + ( includeEikonal ? 2.*beta*lam : 0. )); s4 =-x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta); } } } s1 *= -fact; s2 *= -fact; s3 *= -fact; s4 *= -fact; return SpinorWaveFunction(pout,out,s1,s2,s3,s4); } SpinorBarWaveFunction FFVVertex::evaluateSmall(Energy2 q2,int iopt, tcPDPtr out, const SpinorBarWaveFunction & sbar, const VectorWaveFunction & vec, unsigned int fhel, unsigned int vhel, double ctheta, double phi, double stheta, bool includeEikonal, SmallAngleDirection direction, Energy mass, Energy) { assert(fhel <= 1); assert( vhel == 0 || vhel == 2 ); SpinorBarWaveFunction output; // first calculate the couplings setCoupling(q2,out,sbar.particle(),vec.particle()); Complex ii(0.,1.); if(mass < ZERO) mass = iopt==5 ? ZERO : out->mass(); Lorentz5Momentum pout = sbar.momentum()+vec.momentum(); assert(sbar.direction()!=intermediate); // helicity of the boson double lam = double(vhel)-1.; // energies and velocities Energy Ef = abs(sbar.momentum().e()); Energy Eg = abs(vec .momentum().e()); double x = Eg/Ef; double beta = sqrt(1.-sqr(mass/Ef)); // calculation of propagator accurate as beta->1 and theta -> 0 Energy2 dot = 2.*Ef*Eg*(sqr(mass/Ef)/(1.+beta)*ctheta + sqr(stheta)/(1.+ctheta) ); Complex fact= norm()*(0.5*left()+0.5*right())*UnitRemoval::E2/dot; // calculation of the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); Complex ephig = cos(phi )+ii*sin(phi ); double dm = mass*UnitRemoval::InvE; double dEf = Ef*UnitRemoval::InvE; double rE = sqrt(.5*dEf); // u-type spinor if(sbar.wave().Type()==SpinorType::u) { // fermion along +z if(direction==PostiveZDirection) { if(fhel==0) { s1 = x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta); s2 = rE*dm/sqrt(1.+beta)*stheta* (+x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. )); s3 = -x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta); s4 = rE*dEf*sqrt(1.+beta)*stheta* (+x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. )); } else if(fhel==1) { s1 = -rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. )); s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(lam+1.)*sqr(stheta)/(1.+ctheta); s3 =-rE*dm/sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); s4 = -x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta); s4 = rE*dm*(1.+ctheta)/ephig*x*(1.-lam)/sqrt(1.+beta); } } // fermion aling -z else { if(fhel==0) { s1 = rE*dm/sqrt(1.+beta)*stheta* (+x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. )); s2 = -x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta); s3 = rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. )); s4 = -x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta); } else if(fhel==1) { s1 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s2 = rE*dEf*sqrt(1.+beta)*stheta* (+x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. )); s3 =-x*rE*dm/ephig*(lam-1.)*(1.+ctheta)/sqrt(1.+beta); s4 = rE*dm/sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); } } } // v-type spinor else if(sbar.wave().Type()==SpinorType::v) { // anti fermion along +z if(direction==PostiveZDirection) { if(fhel==0) { s1 = -rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. )); s2 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s3 = rE*dm/sqrt(1.+beta)*stheta* ( x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. )); s4 =-x*rE*dm/ephig*(1.-lam)*(1.+ctheta)/sqrt(1.+beta); } else if(fhel==1) { s1 =-x*rE*dm*ephig*(1.+lam)*(1.+ctheta)/sqrt(1.+beta); s2 =-rE*dm/sqrt(1.+beta)*stheta* ( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. )); s3 =-x*rE*dEf*sqrt(1.+beta)*ephig*(1.-lam)*sqr(stheta)/(1.+ctheta); s4 = rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); } } // anti fermion aling -z else { if(fhel==0) { s1 = -x*rE*dEf*sqrt(1.+beta)/ephig*(1.+lam)*sqr(stheta)/(1.+ctheta); s2 = rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.+lam) - (includeEikonal ? 2.*beta*lam : 0. )); s3 = x*rE*dm/ephig*(-1.+lam)*(1.+ctheta)/sqrt(1.+beta); s4 =-rE*dm/sqrt(1.+beta)*stheta* (+x*(1.-lam) - (includeEikonal ? 2.*beta*lam : 0. )); } else if(fhel==1) { s1 =-rE*dm/sqrt(1.+beta)*stheta* ( x*(1.+lam) + (includeEikonal ? 2.*beta*lam : 0. )); s2 = x*rE*dm*ephig*(lam+1.)*(1.+ctheta)/sqrt(1.+beta); s3 = rE*dEf*sqrt(1.+beta)*stheta* ( x*(1.-lam) + (includeEikonal ? 2.*beta*lam : 0. )); s4 =-x*rE*dEf*sqrt(1.+beta)*ephig*(lam-1.)*sqr(stheta)/(1.+ctheta); } } } s1 *= -fact; s2 *= -fact; s3 *= -fact; s4 *= -fact; return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4); } diff --git a/Helicity/Vertex/Vector/GeneralFFVVertex.cc b/Helicity/Vertex/Vector/GeneralFFVVertex.cc --- a/Helicity/Vertex/Vector/GeneralFFVVertex.cc +++ b/Helicity/Vertex/Vector/GeneralFFVVertex.cc @@ -1,348 +1,348 @@ // -*- C++ -*- // // GeneralFFVVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the GeneralFFVVertex class. // #include "GeneralFFVVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // Definition of the static class description member // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGGeneralFFVVertex("ThePEG::GeneralFFVVertex", "libThePEG.so"); void GeneralFFVVertex::Init() { static ClassDocumentation documentation ("The GeneralFFVVertex class implements the helicity amplitude" "calculations for a fermion-fantifermion gauge boson vertex. Any " "implementation of such a vertex should inherit from in and implement" " the virtual setCoupling member to calculate the coupling"); } // evalulate the full vertex Complex GeneralFFVVertex::evaluate(Energy2 q2, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, const VectorWaveFunction & vec) { // first calculate the couplings setCoupling(q2,sp.particle(),sbar.particle(),vec.particle()); const Complex ii(0.,1.); const complex zero(ZERO,ZERO); // useful combinations of the polarization vector components Complex e0p3=vec.t()+vec.z(); Complex e0m3=vec.t()-vec.z(); Complex e1p2=vec.x()+ii*vec.y(); Complex e1m2=vec.x()-ii*vec.y(); complex p0p3=vec. e()+vec.pz(); complex p0m3=vec. e()-vec.pz(); complex p1p2=vec.px()+ii*vec.py(); complex p1m2=vec.px()-ii*vec.py(); // get both the spinors in the same representation (using the default choice) Complex vertex(0.); // first the left piece as this is virtually always needed if(_left!=0.) { vertex = _left*(+sbar.s3()*(sp.s1()*e0p3+sp.s2()*e1m2) +sbar.s4()*(sp.s1()*e1p2+sp.s2()*e0m3)); } // then the right piece (often not needed eg W vertex) if(_right!=0.) { vertex += _right*(+sbar.s1()*(sp.s3()*e0m3-sp.s4()*e1m2) -sbar.s2()*(sp.s3()*e1p2-sp.s4()*e0p3)); } // left sigma piece if(_leftSigma!=zero) { - vertex += -ii * _leftSigma * - (sbar.s1()*sp.s1()*(-vec. e()*vec.z() + vec.pz()*vec.t() - -ii*(vec.py()*vec.x()-vec.px()*vec.y()))+ - sbar.s2()*sp.s1()*( p1p2*e0p3 - p0p3*e1p2 ) + - sbar.s1()*sp.s2()*( + p1m2*e0m3 - p0m3*e1m2 ) + - sbar.s2()*sp.s2()*( vec. e()*vec.z() - vec.pz()*vec.t() - -ii*(vec.px()*vec.y()-vec.py()*vec.x()))); + vertex -= Complex(ii * _leftSigma * + (sbar.s1()*sp.s1()*(-vec. e()*vec.z() + vec.pz()*vec.t() + -ii*(vec.py()*vec.x()-vec.px()*vec.y()))+ + sbar.s2()*sp.s1()*( p1p2*e0p3 - p0p3*e1p2 ) + + sbar.s1()*sp.s2()*( + p1m2*e0m3 - p0m3*e1m2 ) + + sbar.s2()*sp.s2()*( vec. e()*vec.z() - vec.pz()*vec.t() + -ii*(vec.px()*vec.y()-vec.py()*vec.x())))); } // right sigma piece if(_rightSigma!=zero) { - vertex += ii * _rightSigma * - (sbar.s3()*sp.s3()*(-vec.e()*vec.z()+vec.pz()*vec.t() - -ii*(vec.px()*vec.y()-vec.py()*vec.x()))+ - sbar.s4()*sp.s3()*( p1p2*e0m3 - p0m3*e1p2 ) + - sbar.s3()*sp.s4()*( p1m2*e0p3 - p0p3*e1m2 ) + - sbar.s4()*sp.s4()*( vec.e()*vec.z()-vec.pz()*vec.t() - +ii*(vec.px()*vec.y()-vec.py()*vec.x()))); + vertex += Complex(ii * _rightSigma * + (sbar.s3()*sp.s3()*(-vec.e()*vec.z()+vec.pz()*vec.t() + -ii*(vec.px()*vec.y()-vec.py()*vec.x()))+ + sbar.s4()*sp.s3()*( p1p2*e0m3 - p0m3*e1p2 ) + + sbar.s3()*sp.s4()*( p1m2*e0p3 - p0p3*e1m2 ) + + sbar.s4()*sp.s4()*( vec.e()*vec.z()-vec.pz()*vec.t() + +ii*(vec.px()*vec.y()-vec.py()*vec.x())))); } vertex *= ii; // final factors return vertex*norm(); } // evaluate an off-shell spinor SpinorWaveFunction GeneralFFVVertex::evaluate(Energy2 q2, int iopt,tcPDPtr out, const SpinorWaveFunction & sp, const VectorWaveFunction &vec, complex mass, complex width) { // extract the pointers to the particle data objects tcPDPtr Psp=sp.particle(); tcPDPtr Pvec=vec.particle(); // first calculate the couplings setCoupling(q2,Psp,out,Pvec); const Complex ii(0.,1.); const complex zero(ZERO,ZERO); // work out the momentum of the off-shell particle Lorentz5Momentum pout = sp.momentum()+vec.momentum(); // now evaluate the contribution // polarization components Complex e0p3 = vec.t() + vec.z(); Complex e0m3 = vec.t() - vec.z(); Complex e1p2 = vec.x()+ii*vec.y(); Complex e1m2 = vec.x()-ii*vec.y(); // overall factor Energy2 p2 = pout.m2(); Complex fact=-normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = iopt==5 ? ZERO : out->mass(); complex p0p3 = pout.e() + pout.z(); complex p0m3 = pout.e() - pout.z(); complex p1p2 = pout.x() + ii*pout.y(); complex p1m2 = pout.x() - ii*pout.y(); // complex nos for for the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); Complex a1 = fact*( sp.s3()*e0m3-sp.s4()*e1m2); Complex a2 = fact*(-sp.s3()*e1p2+sp.s4()*e0p3); Complex a3 = fact*( sp.s1()*e0p3+sp.s2()*e1m2); Complex a4 = fact*( sp.s1()*e1p2+sp.s2()*e0m3); // left piece if(_left!=0.) { - s1 +=UnitRemoval::InvE * _left * (p0m3*a3-p1m2*a4); - s2 +=UnitRemoval::InvE * _left * (-p1p2*a3+p0p3*a4); - s3 +=UnitRemoval::InvE * _left * a3*mass; - s4 +=UnitRemoval::InvE * _left * a4*mass; + s1 += Complex(UnitRemoval::InvE * _left * (p0m3*a3-p1m2*a4)); + s2 += Complex(UnitRemoval::InvE * _left * (-p1p2*a3+p0p3*a4)); + s3 += Complex(UnitRemoval::InvE * _left * a3*mass); + s4 += Complex(UnitRemoval::InvE * _left * a4*mass); } // right piece if(_right!=0.) { - s1 +=UnitRemoval::InvE * _right * a1*mass; - s2 +=UnitRemoval::InvE * _right * a2*mass; - s3 +=UnitRemoval::InvE * _right * (p0p3*a1+p1m2*a2); - s4 +=UnitRemoval::InvE * _right * (p1p2*a1+p0m3*a2); + s1 += Complex(UnitRemoval::InvE * _right * a1*mass); + s2 += Complex(UnitRemoval::InvE * _right * a2*mass); + s3 += Complex(UnitRemoval::InvE * _right * (p0p3*a1+p1m2*a2)); + s4 += Complex(UnitRemoval::InvE * _right * (p1p2*a1+p0m3*a2)); } complex p0p3b = vec. e() + vec.pz(); complex p0m3b = vec. e() - vec.pz(); complex p1p2b = vec.px() + ii*vec.py(); complex p1m2b = vec.px() - ii*vec.py(); complex b1 = fact*( sp.s3()*p0m3b - sp.s4()*p1m2b); complex b2 = fact*(-sp.s3()*p1p2b + sp.s4()*p0p3b); complex b3 = fact*( sp.s1()*p0p3b + sp.s2()*p1m2b); complex b4 = fact*( sp.s1()*p1p2b + sp.s2()*p0m3b); // left sigma piece if(_leftSigma!=zero) { - s1 += -0.5*ii*UnitRemoval::InvE *mass*_leftSigma* - ( - a3*p0m3b + a4*p1m2b + b3*e0m3 - b4*e1m2); - s2 += 0.5*ii*UnitRemoval::InvE *mass*_leftSigma* - ( - a3*p1p2b + a4*p0p3b + b3*e1p2 - b4*e0p3); - s3 += -0.5*ii*UnitRemoval::InvE*_leftSigma* - ( + p0p3*( - a3*p0m3b + a4*p1m2b + b3*e0m3 - b4*e1m2 ) - + p1m2*( a3*p1p2b - a4*p0p3b - b3*e1p2 + b4*e0p3 ) ); - s4 += 0.5*ii*UnitRemoval::InvE*_leftSigma* - ( + p1p2*( + a3*p0m3b - a4*p1m2b - b3*e0m3 + b4*e1m2 ) - + p0m3*( - a3*p1p2b + a4*p0p3b + b3*e1p2 - b4*e0p3 ) ); + s1 += Complex(-0.5*ii*UnitRemoval::InvE *mass*_leftSigma* + ( - a3*p0m3b + a4*p1m2b + b3*e0m3 - b4*e1m2)); + s2 += Complex( 0.5*ii*UnitRemoval::InvE *mass*_leftSigma* + ( - a3*p1p2b + a4*p0p3b + b3*e1p2 - b4*e0p3)); + s3 += Complex(-0.5*ii*UnitRemoval::InvE*_leftSigma* + ( + p0p3*( - a3*p0m3b + a4*p1m2b + b3*e0m3 - b4*e1m2 ) + + p1m2*( a3*p1p2b - a4*p0p3b - b3*e1p2 + b4*e0p3 ) )); + s4 += Complex( 0.5*ii*UnitRemoval::InvE*_leftSigma* + ( + p1p2*( + a3*p0m3b - a4*p1m2b - b3*e0m3 + b4*e1m2 ) + + p0m3*( - a3*p1p2b + a4*p0p3b + b3*e1p2 - b4*e0p3 ) )); } // right sigma piece if(_rightSigma!=zero) { - s1 += 0.5*ii*UnitRemoval::InvE *_rightSigma* - ( + p0m3*( + a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 ) - + p1m2*( - a1*p1p2b - a2*p0m3b + b1*e1p2 + b2*e0m3 ) ); - s2 += -0.5*ii*UnitRemoval::InvE *_rightSigma* - ( + p1p2*( + a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 ) - + p0p3*( - a1*p1p2b - a2*p0m3b + b1*e1p2 + b2*e0m3 ) ); - s3 += 0.5*ii*UnitRemoval::InvE *mass*_rightSigma* - ( a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 ); - s4 += -0.5*ii*UnitRemoval::InvE *mass*_rightSigma* - ( -a1*p1p2b - a2*p0m3b + b2*e0m3 + b1*e1p2 ); + s1 +=Complex( 0.5*ii*UnitRemoval::InvE *_rightSigma* + ( + p0m3*( + a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 ) + + p1m2*( - a1*p1p2b - a2*p0m3b + b1*e1p2 + b2*e0m3 ) )); + s2 +=Complex( -0.5*ii*UnitRemoval::InvE *_rightSigma* + ( + p1p2*( + a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 ) + + p0p3*( - a1*p1p2b - a2*p0m3b + b1*e1p2 + b2*e0m3 ) )); + s3 += Complex( 0.5*ii*UnitRemoval::InvE *mass*_rightSigma* + ( a1*p0p3b + a2*p1m2b - b1*e0p3 - b2*e1m2 )); + s4 +=Complex( -0.5*ii*UnitRemoval::InvE *mass*_rightSigma* + ( -a1*p1p2b - a2*p0m3b + b2*e0m3 + b1*e1p2 )); } // return the wavefunction return SpinorWaveFunction(pout,out,s1,s2,s3,s4); } // evaluate an off-shell SpinorBar SpinorBarWaveFunction GeneralFFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const SpinorBarWaveFunction & sbar, const VectorWaveFunction & vec, complex mass, complex width) { // first calculate the couplings setCoupling(q2,out,sbar.particle(),vec.particle()); const Complex ii(0.,1.); const complex zero(ZERO,ZERO); // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+vec.momentum(); // now evaluate the contribution // polarization components Complex e0p3=vec.t() + vec.z(); Complex e0m3=vec.t() - vec.z(); Complex e1p2=vec.x()+ii*vec.y(); Complex e1m2=vec.x()-ii*vec.y(); // overall factor Energy2 p2 = pout.m2(); Complex fact=-normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); complex p1p2=pout.x() + ii*pout.y(); complex p1m2=pout.x() - ii*pout.y(); complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); // complex numbers for the spinor Complex s1(0.),s2(0.),s3(0.),s4(0.); Complex a1 = fact*( sbar.s3()*e0p3+sbar.s4()*e1p2); Complex a2 = fact*( sbar.s3()*e1m2+sbar.s4()*e0m3); Complex a3 = fact*( sbar.s1()*e0m3-sbar.s2()*e1p2); Complex a4 = fact*(-sbar.s1()*e1m2+sbar.s2()*e0p3); complex p0p3b = vec. e() + vec.pz(); complex p0m3b = vec. e() - vec.pz(); complex p1p2b = vec.px() + ii*vec.py(); complex p1m2b = vec.px() - ii*vec.py(); complex b1 = fact*( sbar.s3()*p0p3b+sbar.s4()*p1p2b); complex b2 = fact*( sbar.s3()*p1m2b+sbar.s4()*p0m3b); complex b3 = fact*( sbar.s1()*p0m3b-sbar.s2()*p1p2b); complex b4 = fact*(-sbar.s1()*p1m2b+sbar.s2()*p0p3b); // left piece if(_left!=0.) { - s1 += UnitRemoval::InvE*_left*a1*mass; - s2 += UnitRemoval::InvE*_left*a2*mass; - s3 += UnitRemoval::InvE*_left*(-p0m3*a1+p1p2*a2); - s4 += UnitRemoval::InvE*_left*(+p1m2*a1-p0p3*a2); + s1 += Complex(UnitRemoval::InvE*_left*a1*mass); + s2 += Complex(UnitRemoval::InvE*_left*a2*mass); + s3 += Complex(UnitRemoval::InvE*_left*(-p0m3*a1+p1p2*a2)); + s4 += Complex(UnitRemoval::InvE*_left*(+p1m2*a1-p0p3*a2)); } // right piece if(_right!=0.) { - s1 += UnitRemoval::InvE*_right*(-p0p3*a3-p1p2*a4); - s2 += UnitRemoval::InvE*_right*(-p1m2*a3-p0m3*a4); - s3 += UnitRemoval::InvE*_right*a3*mass; - s4 += UnitRemoval::InvE*_right*a4*mass; + s1 += Complex(UnitRemoval::InvE*_right*(-p0p3*a3-p1p2*a4)); + s2 += Complex(UnitRemoval::InvE*_right*(-p1m2*a3-p0m3*a4)); + s3 += Complex(UnitRemoval::InvE*_right*a3*mass); + s4 += Complex(UnitRemoval::InvE*_right*a4*mass); } // left sigma piece if(_leftSigma!=zero) { - s1 += 0.5*ii*UnitRemoval::InvE*_leftSigma*mass* - ( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2); - s2 += -0.5*ii*UnitRemoval::InvE*_leftSigma*mass* - ( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3); - s3 += -0.5*ii*UnitRemoval::InvE*_leftSigma* - (+p0m3*( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2) - +p1p2*( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3)); - s4 += +0.5*ii*UnitRemoval::InvE*_leftSigma* - (+p1m2*( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2) - +p0p3*( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3)); + s1 +=Complex( 0.5*ii*UnitRemoval::InvE*_leftSigma*mass* + ( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2)); + s2 +=Complex( -0.5*ii*UnitRemoval::InvE*_leftSigma*mass* + ( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3)); + s3 +=Complex( -0.5*ii*UnitRemoval::InvE*_leftSigma* + (+p0m3*( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2) + +p1p2*( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3))); + s4 +=Complex( +0.5*ii*UnitRemoval::InvE*_leftSigma* + (+p1m2*( + a3*p0p3b + a4*p1p2b - b3*e0p3 - b4*e1p2) + +p0p3*( - a3*p1m2b - a4*p0m3b + b3*e1m2 + b4*e0m3))); } // right sigma piece if(_rightSigma!=zero) { - s1 += +0.5*ii*UnitRemoval::InvE*_rightSigma* - (+p0p3*( - a1*p0m3b + a2*p1p2b + b1*e0m3 - b2*e1p2) - +p1p2*( + a1*p1m2b - a2*p0p3b - b1*e1m2 + b2*e0p3)); - s2 += -0.5*ii*UnitRemoval::InvE*_rightSigma* - (+p1m2*( + a1*p0m3b - a2*p1p2b - b1*e0m3 + b2*e1p2) - +p0m3*( - a1*p1m2b + a2*p0p3b + b1*e1m2 - b2*e0p3)); - s3 += -0.5*ii*UnitRemoval::InvE*_rightSigma*mass* - ( - a1*p0m3b + a2*p1p2b + b1*e0m3 - b2*e1p2 ); - s4 += 0.5*ii*UnitRemoval::InvE*_rightSigma*mass* - ( - a1*p1m2b + a2*p0p3b + b1*e1m2 - b2*e0p3 ); + s1 +=Complex( +0.5*ii*UnitRemoval::InvE*_rightSigma* + (+p0p3*( - a1*p0m3b + a2*p1p2b + b1*e0m3 - b2*e1p2) + +p1p2*( + a1*p1m2b - a2*p0p3b - b1*e1m2 + b2*e0p3))); + s2 +=Complex( -0.5*ii*UnitRemoval::InvE*_rightSigma* + (+p1m2*( + a1*p0m3b - a2*p1p2b - b1*e0m3 + b2*e1p2) + +p0m3*( - a1*p1m2b + a2*p0p3b + b1*e1m2 - b2*e0p3))); + s3 +=Complex( -0.5*ii*UnitRemoval::InvE*_rightSigma*mass* + ( - a1*p0m3b + a2*p1p2b + b1*e0m3 - b2*e1p2 )); + s4 +=Complex( 0.5*ii*UnitRemoval::InvE*_rightSigma*mass* + ( - a1*p1m2b + a2*p0p3b + b1*e1m2 - b2*e0p3 )); } return SpinorBarWaveFunction(pout,out,s1,s2,s3,s4); } // off-shell vector VectorWaveFunction GeneralFFVVertex::evaluate(Energy2 q2,int iopt,tcPDPtr out, const SpinorWaveFunction & sp, const SpinorBarWaveFunction & sbar, complex mass, complex width) { // first calculate the couplings setCoupling(q2,sp.particle(),sbar.particle(),out); const Complex ii(0.,1.); const complex zero(ZERO,ZERO); // work out the momentum of the off-shell particle Lorentz5Momentum pout = sbar.momentum()+sp.momentum(); // overall factor Energy2 p2 = pout.m2(); Complex fact = normPropagator(iopt,p2,out,mass,width); // momentum components if(mass.real() < ZERO) mass = (iopt==5) ? ZERO : out->mass(); complex mass2 = sqr(mass); // the vector for the fermion-antifermion Complex vec[4]; complex p1p2=pout.x() + ii*pout.y(); complex p1m2=pout.x() - ii*pout.y(); complex p0p3=pout.e() + pout.z(); complex p0m3=pout.e() - pout.z(); // left coupling if(_left!=0.) { vec[0] = -_left*(sbar.s3()*sp.s2()+sbar.s4()*sp.s1()); vec[1] = ii*_left*(sbar.s3()*sp.s2()-sbar.s4()*sp.s1()); vec[2] = -_left*(sbar.s3()*sp.s1()-sbar.s4()*sp.s2()); vec[3] = _left*(sbar.s3()*sp.s1()+sbar.s4()*sp.s2()); } // right coupling if(_right!=0.) { vec[0] += +_right*(sbar.s1()*sp.s4()+sbar.s2()*sp.s3()); vec[1] += -ii*_right*(sbar.s1()*sp.s4()-sbar.s2()*sp.s3()); vec[2] += +_right*(sbar.s1()*sp.s3()-sbar.s2()*sp.s4()); vec[3] += +_right*(sbar.s1()*sp.s3()+sbar.s2()*sp.s4()); } // left sigma if(_leftSigma==zero) { - vec[0] += -0.5*ii*_leftSigma* - (+sp.s1()*sbar.s1()*(p1p2-p1m2)+2.*sp.s1()*sbar.s2()*p0p3 - -sp.s2()*sbar.s2()*(p1p2-p1m2)+2.*sp.s2()*sbar.s1()*p0m3); - vec[1] += -0.5 *_leftSigma* - (+sp.s1()*sbar.s1()*(p1m2-p1p2)+2.*sp.s1()*sbar.s2()*p0p3 - +sp.s2()*sbar.s2()*(p1p2+p1m2)-2.*sp.s2()*sbar.s1()*p0m3); - vec[2] += -0.5*ii*_leftSigma* - (+sp.s1()*sbar.s1()*(p0p3+p0m3)-2.*sp.s1()*sbar.s2()*p1p2 - -sp.s2()*sbar.s2()*(p0p3+p0m3)+2.*sp.s2()*sbar.s1()*p1m2); - vec[3] += 0.5*ii*_leftSigma* - (-sp.s1()*sbar.s1()*(p0p3-p0m3)-2.*sp.s1()*sbar.s2()*p1p2 - +sp.s2()*sbar.s2()*(p0p3-p0m3)-2.*sp.s2()*sbar.s1()*p1m2); + vec[0] +=Complex( -0.5*ii*_leftSigma* + (+sp.s1()*sbar.s1()*(p1p2-p1m2)+2.*sp.s1()*sbar.s2()*p0p3 + -sp.s2()*sbar.s2()*(p1p2-p1m2)+2.*sp.s2()*sbar.s1()*p0m3)); + vec[1] +=Complex( -0.5 *_leftSigma* + (+sp.s1()*sbar.s1()*(p1m2-p1p2)+2.*sp.s1()*sbar.s2()*p0p3 + +sp.s2()*sbar.s2()*(p1p2+p1m2)-2.*sp.s2()*sbar.s1()*p0m3)); + vec[2] +=Complex( -0.5*ii*_leftSigma* + (+sp.s1()*sbar.s1()*(p0p3+p0m3)-2.*sp.s1()*sbar.s2()*p1p2 + -sp.s2()*sbar.s2()*(p0p3+p0m3)+2.*sp.s2()*sbar.s1()*p1m2)); + vec[3] +=Complex( 0.5*ii*_leftSigma* + (-sp.s1()*sbar.s1()*(p0p3-p0m3)-2.*sp.s1()*sbar.s2()*p1p2 + +sp.s2()*sbar.s2()*(p0p3-p0m3)-2.*sp.s2()*sbar.s1()*p1m2)); } // right sigma if(_rightSigma==zero) { - vec[0] += 0.5*ii*_rightSigma* - (-sp.s3()*sbar.s3()*(p1p2-p1m2)+2.*sp.s3()*sbar.s4()*p0m3 - +sp.s4()*sbar.s4()*(p1p2-p1m2)+2.*sp.s4()*sbar.s3()*p0p3); - vec[1] += 0.5* _rightSigma* - (-sp.s3()*sbar.s3()*(p1p2+p1m2)-2.*sp.s3()*sbar.s4()*p0m3 - +sp.s4()*sbar.s4()*(p1p2+p1m2)+2.*sp.s4()*sbar.s3()*p0p3); - vec[2] += 0.5*ii*_rightSigma* - (+sp.s3()*sbar.s3()*(p0p3+p0m3)+2.*sp.s3()*sbar.s4()*p1p2 - -sp.s4()*sbar.s4()*(p0p3+p0m3)-2.*sp.s4()*sbar.s3()*p1m2); - vec[3] += -0.5*ii*_rightSigma* - (-sp.s3()*sbar.s3()*(p0p3-p0m3)-2.*sp.s3()*sbar.s4()*p1p2 - +sp.s4()*sbar.s4()*(p0p3-p0m3)-2.*sp.s4()*sbar.s3()*p1m2); + vec[0] +=Complex( 0.5*ii*_rightSigma* + (-sp.s3()*sbar.s3()*(p1p2-p1m2)+2.*sp.s3()*sbar.s4()*p0m3 + +sp.s4()*sbar.s4()*(p1p2-p1m2)+2.*sp.s4()*sbar.s3()*p0p3)); + vec[1] +=Complex( 0.5* _rightSigma* + (-sp.s3()*sbar.s3()*(p1p2+p1m2)-2.*sp.s3()*sbar.s4()*p0m3 + +sp.s4()*sbar.s4()*(p1p2+p1m2)+2.*sp.s4()*sbar.s3()*p0p3)); + vec[2] +=Complex( 0.5*ii*_rightSigma* + (+sp.s3()*sbar.s3()*(p0p3+p0m3)+2.*sp.s3()*sbar.s4()*p1p2 + -sp.s4()*sbar.s4()*(p0p3+p0m3)-2.*sp.s4()*sbar.s3()*p1m2)); + vec[3] +=Complex( -0.5*ii*_rightSigma* + (-sp.s3()*sbar.s3()*(p0p3-p0m3)-2.*sp.s3()*sbar.s4()*p1p2 + +sp.s4()*sbar.s4()*(p0p3-p0m3)-2.*sp.s4()*sbar.s3()*p1m2)); } // massless boson if(mass.real()==ZERO) { for(int ix=0;ix<4;++ix){vec[ix]*=fact;} } // massive boson else { complex dot = ( pout.e() *vec[3] -pout.x()*vec[0] -pout.y()*vec[1] -pout.z()*vec[2])/mass2; vec[0]=fact*(vec[0]-dot*pout.x()); vec[1]=fact*(vec[1]-dot*pout.y()); vec[2]=fact*(vec[2]-dot*pout.z()); vec[3]=fact*(vec[3]-dot*pout.e()); } return VectorWaveFunction(pout,out,vec[0],vec[1],vec[2],vec[3]); } diff --git a/Helicity/Vertex/Vector/VVVVVertex.cc b/Helicity/Vertex/Vector/VVVVVertex.cc --- a/Helicity/Vertex/Vector/VVVVVertex.cc +++ b/Helicity/Vertex/Vector/VVVVVertex.cc @@ -1,313 +1,313 @@ // -*- C++ -*- // // VVVVVertex.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the VVVVVertex class. // #include "VVVVVertex.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" using namespace ThePEG; using namespace Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeThePEGVVVVVertex("ThePEG::VVVVVertex", "libThePEG.so"); void VVVVVertex::Init() { static ClassDocumentation documentation ("The VVVVVertex class is the implementation of the 4-vector vertex"); } // calculate the vertex Complex VVVVVertex::evaluate(Energy2 q2 , int iopt, const VectorWaveFunction & vec1, const VectorWaveFunction & vec2, const VectorWaveFunction & vec3, const VectorWaveFunction & vec4) { // workout the coupling setCoupling(q2,vec1.particle(),vec2.particle(), vec3.particle(),vec4.particle()); Complex vertex,ii(0.,1.); // calculate the vertex // QCD type vertex assert(_itype>=1&&_itype<=2); if(_itype==1) { // dot products we need Complex dotv1v2 = vec1.wave().dot(vec2.wave()); Complex dotv3v4 = vec3.wave().dot(vec4.wave()); Complex dotv1v4 = vec1.wave().dot(vec4.wave()); Complex dotv2v3 = vec3.wave().dot(vec2.wave()); // first the 4-point part of the vertex vertex = dotv1v2*dotv3v4-dotv1v4*dotv2v3; // now the virtual gluon exchange if needed if(iopt!=0) { // dot products Complex dotv1v3 = vec1.wave().dot(vec3.wave()); Complex dotv2v4 = vec4.wave().dot(vec2.wave()); complex dotv1p13 = vec1.wave().dot(2.*vec3.momentum() + vec1.momentum()); complex dotv2p24 = vec2.wave().dot(2.*vec4.momentum() + vec2.momentum()); complex dotv3p13 = vec3.wave().dot(2.*vec1.momentum() + vec3.momentum()); complex dotv4p24 = vec4.wave().dot(2.*vec2.momentum() + vec4.momentum()); LorentzPolarizationVectorE veca = dotv3p13*vec1.wave() - dotv1p13*vec3.wave() + dotv1v3*(vec3.momentum()-vec1.momentum()); LorentzPolarizationVectorE vecb = dotv4p24*vec2.wave() - dotv2p24*vec4.wave() + dotv2v4*(vec4.momentum()-vec2.momentum()); InvEnergy2 numerator = 1./(vec1.momentum()+vec3.momentum()).m2(); - vertex += numerator*veca.dot(vecb); + vertex += Complex(numerator*veca.dot(vecb)); } } // EW type vertex else if(_itype==2) { Complex dotv1v2 = vec1.wave().dot(vec2.wave()); Complex dotv1v3 = vec1.wave().dot(vec3.wave()); Complex dotv1v4 = vec1.wave().dot(vec4.wave()); Complex dotv2v3 = vec2.wave().dot(vec3.wave()); Complex dotv2v4 = vec2.wave().dot(vec4.wave()); Complex dotv3v4 = vec3.wave().dot(vec4.wave()); // evaluate the vertex // need to sort the order out here if(( _iorder[0]==0 && _iorder[1]==1 && _iorder[2]==2 && _iorder[3]==3)|| ( _iorder[0]==1 && _iorder[1]==0 && _iorder[2]==2 && _iorder[3]==3)|| ( _iorder[0]==0 && _iorder[1]==1 && _iorder[2]==3 && _iorder[3]==2)|| ( _iorder[0]==1 && _iorder[1]==0 && _iorder[2]==3 && _iorder[3]==2)|| ( _iorder[0]==2 && _iorder[1]==3 && _iorder[2]==0 && _iorder[3]==1)|| ( _iorder[0]==2 && _iorder[1]==3 && _iorder[2]==1 && _iorder[3]==0)|| ( _iorder[0]==3 && _iorder[1]==2 && _iorder[2]==0 && _iorder[3]==1)|| ( _iorder[0]==3 && _iorder[1]==2 && _iorder[2]==1 && _iorder[3]==0)) { // contact term vertex = 2.*dotv1v2*dotv3v4-dotv1v3*dotv2v4-dotv1v4*dotv2v3; // now for the u- and t-channel terms if needed if(iopt!=0) { // dot products of momenta and wavefunction complex dotv1p13 = +vec1.wave().dot(vec1.momentum() + 2.*vec3.momentum() ); complex dotv1p14 = +vec1.wave().dot(vec1.momentum() + 2.*vec4.momentum() ); complex dotv2p23 = +vec2.wave().dot(vec2.momentum() + 2.*vec3.momentum() ); complex dotv2p24 = +vec2.wave().dot(vec2.momentum() + 2.*vec4.momentum() ); complex dotv3p31 = +vec3.wave().dot(vec3.momentum() + 2.*vec1.momentum() ); complex dotv3p32 = +vec3.wave().dot(vec3.momentum() + 2.*vec2.momentum() ); complex dotv4p41 = +vec4.wave().dot(vec4.momentum() + 2.*vec1.momentum() ); complex dotv4p42 = +vec4.wave().dot(vec4.momentum() + 2.*vec2.momentum() ); LorentzPolarizationVectorE ja = (vec3.momentum()-vec1.momentum())*dotv1v3 +dotv3p31*vec1.wave()-dotv1p13*vec3.wave(); LorentzPolarizationVectorE jb = (vec4.momentum()-vec2.momentum())*dotv2v4 +dotv4p42*vec2.wave()-dotv2p24*vec4.wave(); LorentzPolarizationVectorE jc = (vec4.momentum()-vec1.momentum())*dotv1v4 +dotv4p41*vec1.wave()-dotv1p14*vec4.wave(); LorentzPolarizationVectorE jd = (vec3.momentum()-vec2.momentum())*dotv2v3 +dotv3p32*vec2.wave()-dotv2p23*vec3.wave(); // dot products of these vectors complex dotjajb = ja.dot(jb); complex dotjcjd = jc.dot(jd); complex dotjaq = ja.dot(vec1.momentum()+vec3.momentum()); complex dotjbq = jb.dot(vec1.momentum()+vec3.momentum()); complex dotjck = jc.dot(vec1.momentum()+vec4.momentum()); complex dotjdk = jd.dot(vec1.momentum()+vec4.momentum()); Energy2 q2 = (vec1.momentum()+vec3.momentum()).m2(); Energy2 k2 = (vec1.momentum()+vec4.momentum()).m2(); // compute the term we need Energy2 mass2; for(int ix=0;ix<2;++ix) { if(_inter[ix]) { mass2 = sqr(_inter[ix]->mass()); if(mass2!=Energy2()) { - vertex += UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,q2,_inter[ix])* - (dotjajb-dotjaq*dotjbq/mass2); - vertex += UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,k2,_inter[ix])* - (dotjcjd-dotjck*dotjdk/mass2); + vertex += Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,q2,_inter[ix])* + (dotjajb-dotjaq*dotjbq/mass2)); + vertex += Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,k2,_inter[ix])* + (dotjcjd-dotjck*dotjdk/mass2)); } else { - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,q2,_inter[ix])*dotjajb; - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,k2,_inter[ix])*dotjcjd; + vertex+= Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,q2,_inter[ix])*dotjajb); + vertex+= Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,k2,_inter[ix])*dotjcjd); } } } } } else if(( _iorder[0]==0 && _iorder[1]==2 && _iorder[2]==1 && _iorder[3]==3)|| ( _iorder[0]==2 && _iorder[1]==0 && _iorder[2]==1 && _iorder[3]==3)|| ( _iorder[0]==0 && _iorder[1]==2 && _iorder[2]==3 && _iorder[3]==1)|| ( _iorder[0]==2 && _iorder[1]==0 && _iorder[2]==3 && _iorder[3]==1)|| ( _iorder[0]==1 && _iorder[1]==3 && _iorder[2]==0 && _iorder[3]==2)|| ( _iorder[0]==1 && _iorder[1]==3 && _iorder[2]==2 && _iorder[3]==0)|| ( _iorder[0]==3 && _iorder[1]==1 && _iorder[2]==0 && _iorder[3]==2)|| ( _iorder[0]==3 && _iorder[1]==1 && _iorder[2]==2 && _iorder[3]==0)) { // contact term vertex = 2.*dotv1v3*dotv2v4-dotv1v2*dotv3v4-dotv1v4*dotv2v3; // now for the u- and t-channel terms if needed if(iopt!=0) { // dot products of momenta and wavefunction complex dotv1p12 = vec1.wave().dot(vec1.momentum() + 2.*vec2.momentum() ); complex dotv1p14 = vec1.wave().dot(vec1.momentum() + 2.*vec4.momentum() ); complex dotv3p32 = vec3.wave().dot(vec3.momentum() + 2.*vec2.momentum() ); complex dotv3p34 = vec3.wave().dot(vec3.momentum() + 2.*vec4.momentum() ); complex dotv2p21 = vec2.wave().dot(vec2.momentum() + 2.*vec1.momentum() ); complex dotv2p23 = vec2.wave().dot(vec2.momentum() + 2.*vec3.momentum() ); complex dotv4p41 = vec4.wave().dot(vec4.momentum() + 2.*vec1.momentum() ); complex dotv4p43 = vec4.wave().dot(vec4.momentum() + 2.*vec3.momentum() ); LorentzPolarizationVectorE ja = (vec2.momentum() - vec1.momentum() )*dotv1v2 + dotv2p21*vec1.wave() - dotv1p12*vec2.wave(); LorentzPolarizationVectorE jb = (vec4.momentum() - vec3.momentum())*dotv3v4 + dotv4p43*vec3.wave() - dotv3p34*vec4.wave(); LorentzPolarizationVectorE jc = (vec4.momentum() - vec1.momentum())*dotv1v4 + dotv4p41*vec1.wave() - dotv1p14*vec4.wave(); LorentzPolarizationVectorE jd = (vec2.momentum() - vec3.momentum())*dotv2v3 + dotv2p23*vec3.wave() - dotv3p32*vec2.wave(); // dot products of these vectors complex dotjajb = ja.dot(jb); complex dotjcjd = jc.dot(jd); complex dotjaq = ja.dot(vec1.momentum()+vec2.momentum()); complex dotjbq = jb.dot(vec1.momentum()+vec2.momentum()); complex dotjck = jc.dot(vec1.momentum()+vec4.momentum()); complex dotjdk = jd.dot(vec1.momentum()+vec4.momentum()); Energy2 q2 = (vec1.momentum()+vec2.momentum()).m2(); Energy2 k2 = (vec1.momentum()+vec4.momentum()).m2(); // compute the term we need Energy2 mass2; for(int ix=0;ix<2;++ix) { if(_inter[ix]) { mass2 = (_inter[ix]->mass())*(_inter[ix]->mass()); if(mass2!=Energy2()) { - vertex+=UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,q2,_inter[ix])* - (dotjajb-dotjaq*dotjbq/mass2); - vertex+=UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,k2,_inter[ix])* - (dotjcjd-dotjck*dotjdk/mass2); + vertex+=Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,q2,_inter[ix])* + (dotjajb-dotjaq*dotjbq/mass2)); + vertex+=Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,k2,_inter[ix])* + (dotjcjd-dotjck*dotjdk/mass2)); } else { - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,q2,_inter[ix])*dotjajb; - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,k2,_inter[ix])*dotjcjd; + vertex+=Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,q2,_inter[ix])*dotjajb); + vertex+=Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,k2,_inter[ix])*dotjcjd); } } } } } else if(( _iorder[0]==0 && _iorder[1]==3 && _iorder[2]==1 && _iorder[3]==2)|| ( _iorder[0]==3 && _iorder[1]==0 && _iorder[2]==1 && _iorder[3]==2)|| ( _iorder[0]==0 && _iorder[1]==3 && _iorder[2]==2 && _iorder[3]==1)|| ( _iorder[0]==3 && _iorder[1]==0 && _iorder[2]==2 && _iorder[3]==1)|| ( _iorder[0]==1 && _iorder[1]==2 && _iorder[2]==0 && _iorder[3]==3)|| ( _iorder[0]==2 && _iorder[1]==1 && _iorder[2]==0 && _iorder[3]==3)|| ( _iorder[0]==1 && _iorder[1]==2 && _iorder[2]==3 && _iorder[3]==0)|| ( _iorder[0]==2 && _iorder[1]==1 && _iorder[2]==3 && _iorder[3]==0)) { // contact term vertex = 2.*dotv1v4*dotv2v3-dotv1v3*dotv2v4-dotv1v2*dotv3v4; // now for the u- and t-channel terms if needed if(iopt!=0) { // dot products of momenta and wavefunction complex dotv1p12 = vec1.wave().dot(vec1.momentum() + 2.*vec2.momentum() ); complex dotv1p13 = vec1.wave().dot(vec1.momentum() + 2.*vec3.momentum() ); complex dotv2p24 = vec2.wave().dot(vec2.momentum() + 2.*vec4.momentum() ); complex dotv2p21 = vec2.wave().dot(vec2.momentum() + 2.*vec1.momentum() ); complex dotv3p31 = vec3.wave().dot(vec3.momentum() + 2.*vec1.momentum() ); complex dotv3p34 = vec3.wave().dot(vec3.momentum() + 2.*vec4.momentum() ); complex dotv4p43 = vec4.wave().dot(vec4.momentum() + 2.*vec3.momentum() ); complex dotv4p42 = vec4.wave().dot(vec4.momentum() + 2.*vec2.momentum() ); LorentzPolarizationVectorE ja = (vec2.momentum()-vec1.momentum())*dotv1v2 +dotv2p21*vec1.wave()-dotv1p12*vec2.wave(); LorentzPolarizationVectorE jb = (vec3.momentum()-vec4.momentum())*dotv3v4 +dotv3p34*vec4.wave()-dotv4p43*vec3.wave(); LorentzPolarizationVectorE jc = (vec3.momentum()-vec1.momentum())*dotv1v3 +dotv3p31*vec1.wave()-dotv1p13*vec3.wave(); LorentzPolarizationVectorE jd = (vec2.momentum()-vec4.momentum())*dotv2v4 +dotv2p24*vec4.wave()-dotv4p42*vec2.wave(); // dot products of these vectors complex dotjajb = ja.dot(jb); complex dotjcjd = jc.dot(jd); complex dotjaq = ja.dot(vec1.momentum()+vec2.momentum()); complex dotjbq = jb.dot(vec1.momentum()+vec2.momentum()); complex dotjck = jc.dot(vec1.momentum()+vec3.momentum()); complex dotjdk = jd.dot(vec1.momentum()+vec3.momentum()); Energy2 q2 = (vec1.momentum()+vec2.momentum()).m2(); Energy2 k2 = (vec1.momentum()+vec3.momentum()).m2(); // compute the term we need Energy2 mass2; for(int ix=0;ix<2;++ix) { if(_inter[ix]) { mass2 = sqr(_inter[ix]->mass()); if(mass2!=Energy2()) { - vertex+=UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,q2,_inter[ix])* - (dotjajb-dotjaq*dotjbq/mass2); - vertex+=UnitRemoval::InvE2 * - _coup[ix]*propagator(iopt,k2,_inter[ix])* - (dotjcjd-dotjck*dotjdk/mass2); + vertex+=Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,q2,_inter[ix])* + (dotjajb-dotjaq*dotjbq/mass2)); + vertex+=Complex(UnitRemoval::InvE2 * + _coup[ix]*propagator(iopt,k2,_inter[ix])* + (dotjcjd-dotjck*dotjdk/mass2)); } else { - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,q2,_inter[ix])*dotjajb; - vertex+=UnitRemoval::InvE2 *_coup[ix]* - propagator(iopt,k2,_inter[ix])*dotjcjd; + vertex+=Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,q2,_inter[ix])*dotjajb); + vertex+=Complex(UnitRemoval::InvE2 *_coup[ix]* + propagator(iopt,k2,_inter[ix])*dotjcjd); } } } } } else throw HelicityConsistencyError() << "Unknown order of particles in " << "VVVVVertex::evaluate()" << Exception::runerror; } // return the answer return -ii*norm()*vertex; } diff --git a/Helicity/WaveFunction/RSSpinorBarWaveFunction.h b/Helicity/WaveFunction/RSSpinorBarWaveFunction.h --- a/Helicity/WaveFunction/RSSpinorBarWaveFunction.h +++ b/Helicity/WaveFunction/RSSpinorBarWaveFunction.h @@ -1,361 +1,361 @@ // -*- C++ -*- // // RSSpinorBarWaveFunction.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_RSSpinorBarWaveFunction_H #define ThePEG_RSSpinorBarWaveFunction_H // // This is the declaration of the RSSpinorBarWaveFunction class. // #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** \ingroup Helicity * * The RSSpinorBarWaveFunction class is designed to * store the wavefunction of a spin-\f$\frac32\f$ particle in a form * suitable for use in helicity amplitude calculations of the matrix * element using a similar philosophy to the FORTRAN HELAS code. * * In addition to storing the barred spinor using the * LorentzRSSpinorBar class it inherits from the * WaveFunctionBase class to provide storage of the * momentum and ParticleData for the fermion. * * This class also contains the code which does the actually * calculation of the barred spinor for an external particle * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$\bar{v}\f$ spinor. * - outgoing calculates a \f$\bar{u}\f$ spinor. * * The barred spinors are calculated using a Clebsch-Gordon decomposition * in the rest-frame for a massive particle and boosted to the lab-frame. * For massless particles the calculation is performed in the lab-frame * (N.B. there are only two helicities \f$\pm\frac32\f$ in this case.) * * N.B. In our convention 0 is the \f$-\frac32\f$ helicity state, * 1 is the \f$-\frac12\f$ helicity state, * 2 is the \f$+\frac12\f$ helicity state * 3 is the \f$+\frac32\f$ helicity state and * * @see WaveFunctionBase * @see LorentzRSSpinorBar * @see HelicityDefinitions * */ class RSSpinorBarWaveFunction: public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor * @param p The momentum. * @param part The ParticleData pointer. * @param xs1 The first spinor component of the \f$x\f$ vector. * @param xs2 The second spinor component of the \f$x\f$ vector. * @param xs3 The third spinor component of the \f$x\f$ vector. * @param xs4 The fourth spinor component of the \f$x\f$ vector. * @param ys1 The first spinor component of the \f$y\f$ vector. * @param ys2 The second spinor component of the \f$y\f$ vector. * @param ys3 The third spinor component of the \f$y\f$ vector. * @param ys4 The fourth spinor component of the \f$y\f$ vector. * @param zs1 The first spinor component of the \f$z\f$ vector. * @param zs2 The second spinor component of the \f$z\f$ vector. * @param zs3 The third spinor component of the \f$z\f$ vector. * @param zs4 The fourth spinor component of the \f$z\f$ vector. * @param ts1 The first spinor component of the \f$t\f$ vector. * @param ts2 The second spinor component of the \f$t\f$ vector. * @param ts3 The third spinor component of the \f$t\f$ vector. * @param ts4 The fourth spinor component of the \f$t\f$ vector. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex xs1, complex xs2, complex xs3, complex xs4, complex ys1, complex ys2, complex ys3, complex ys4, complex zs1, complex zs2, complex zs3, complex zs4, complex ts1, complex ts2, complex ts3, complex ts4) : WaveFunctionBase(p,part), _wf(xs1,xs2,xs3,xs4, ys1,ys2,ys3,ys4, zs1,zs2,zs3,zs4, ts1,ts2,ts3,ts4) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, LorentzRSSpinorBar & wave) : WaveFunctionBase(p,part), _wf(wave) { assert(iSpin()==4); } /** * Constructor, set the particle and the wavefunction. * @param p Particle * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorBarWaveFunction(const tPPtr & p, const LorentzRSSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==4); for (unsigned int i=0; i<4; ++i) for(unsigned int j=0; j<4; ++j) - _wf(i,j)=wave(i,j)*UnitRemoval::InvSqrtE; + _wf(i,j)=Complex(wave(i,j)*UnitRemoval::InvSqrtE); } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1,2,3 as described above.) * @param dir The direction. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel,Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==4); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ RSSpinorBarWaveFunction(Lorentz5Momentum p,tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==4); } /** * Default constructor */ RSSpinorBarWaveFunction() : WaveFunctionBase(), _wf() {} //@} /** * Access to the wavefunction and its components. */ //@{ /** * subscript operator for the wavefunction * Set components by index. */ complex operator ()(int i, int j) const { assert( i>=0 && i<=3 && j>=0 && j<=3 ); return _wf(i,j); } /** * return wavefunction as LorentzRSSpinorBar */ const LorentzRSSpinorBar & wave() const {return _wf;} /** * Get first spinor component for the x vector */ complex xs1() const {return _wf.xs1();} /** * Get second spinor component for the x vector */ complex xs2() const {return _wf.xs2();} /** * Get third spinor component for the x vector */ complex xs3() const {return _wf.xs3();} /** * Get fourth spinor component for the x vector */ complex xs4() const {return _wf.xs4();} /** * Get first spinor component for the y vector */ complex ys1() const {return _wf.ys1();} /** * Get second spinor component for the y vector */ complex ys2() const {return _wf.ys2();} /** * Get third spinor component for the y vector */ complex ys3() const {return _wf.ys3();} /** * Get fourth spinor component for the y vector */ complex ys4() const {return _wf.ys4();} /** * Get first spinor component for the z vector */ complex zs1() const {return _wf.zs1();} /** * Get second spinor component for the z vector */ complex zs2() const {return _wf.zs2();} /** * Get third spinor component for the z vector */ complex zs3() const {return _wf.zs3();} /** * Get fourth spinor component for the z vector */ complex zs4() const {return _wf.zs4();} /** * Get first spinor component for the t vector */ complex ts1() const {return _wf.ts1();} /** * Get second spinor component for the t vector */ complex ts2() const {return _wf.ts2();} /** * Get third spinor component for the t vector */ complex ts3() const {return _wf.ts3();} /** * Get fourth spinor component for the t vector */ complex ts4() const {return _wf.ts4();} //@} /** * reset functions */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1,2,3 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1,2,3 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * storage of the Lorentz RSSpinorBar */ LorentzRSSpinorBar _wf; /// Return wavefunction as LorentzRSSpinorBar LorentzRSSpinorBar dimensionedWf() const { LorentzRSSpinorBar temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) for (unsigned int j=0; j<4; ++j) temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; return temp; } }; } } #endif /* ThePEG_RSSpinorBarWaveFunction_H */ diff --git a/Helicity/WaveFunction/RSSpinorWaveFunction.h b/Helicity/WaveFunction/RSSpinorWaveFunction.h --- a/Helicity/WaveFunction/RSSpinorWaveFunction.h +++ b/Helicity/WaveFunction/RSSpinorWaveFunction.h @@ -1,360 +1,360 @@ // -*- C++ -*- // // RSSpinorWaveFunction.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_RSSpinorWaveFunction_H #define ThePEG_RSSpinorWaveFunction_H // This is the declaration of the RSSpinorWaveFunction class. #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** \ingroup Helicity * * The RSSpinorWaveFunction class is designed to store the wavefunction * of a spin-3/2 particle in a form suitable for use in helicity amplitude * calculations of the matrix element using a similar philosophy to the * FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzRSSpinor class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and ParticleData for the fermion. * * This class also contains the code which does the actually calculation of the * spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$u\f$ spinor. * - outgoing calculates a \f$v\f$ spinor. * * The spinors are calculated using a Clebsch-Gordon decomposition in the rest-frame * for a massive particle and boosted to the lab-frame. For massless particles the * calculation is performed in the lab-frame (N.B. there are only two helicities * \f$\pm\frac32\f$ in this case.) * * N.B. In our convention 0 is the \f$-\frac32\f$ helicity state, * 1 is the \f$-\frac12\f$ helicity state, * 2 is the \f$+\frac12\f$ helicity state * 3 is the \f$+\frac32\f$ helicity state and * * @see WaveFunctionBase * @see LorentzRSSpinor * @see HelicityDefinitions * */ class RSSpinorWaveFunction: public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param xs1 The first spinor component of the \f$x\f$ vector. * @param xs2 The second spinor component of the \f$x\f$ vector. * @param xs3 The third spinor component of the \f$x\f$ vector. * @param xs4 The fourth spinor component of the \f$x\f$ vector. * @param ys1 The first spinor component of the \f$y\f$ vector. * @param ys2 The second spinor component of the \f$y\f$ vector. * @param ys3 The third spinor component of the \f$y\f$ vector. * @param ys4 The fourth spinor component of the \f$y\f$ vector. * @param zs1 The first spinor component of the \f$z\f$ vector. * @param zs2 The second spinor component of the \f$z\f$ vector. * @param zs3 The third spinor component of the \f$z\f$ vector. * @param zs4 The fourth spinor component of the \f$z\f$ vector. * @param ts1 The first spinor component of the \f$t\f$ vector. * @param ts2 The second spinor component of the \f$t\f$ vector. * @param ts3 The third spinor component of the \f$t\f$ vector. * @param ts4 The fourth spinor component of the \f$t\f$ vector. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex xs1, complex xs2, complex xs3, complex xs4, complex ys1, complex ys2, complex ys3, complex ys4, complex zs1, complex zs2, complex zs3, complex zs4, complex ts1, complex ts2, complex ts3, complex ts4) : WaveFunctionBase(p,part), _wf(xs1,xs2,xs3,xs4, ys1,ys2,ys3,ys4, zs1,zs2,zs3,zs4, ts1,ts2,ts3,ts4) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part, const LorentzRSSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The Particle pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorWaveFunction(const tPPtr & p, const LorentzRSSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==4); for (unsigned int i=0; i<4; ++i) for(unsigned int j=0; j<4; ++j) - _wf(i,j)=wave(i,j)*UnitRemoval::InvSqrtE; + _wf(i,j)=Complex(wave(i,j)*UnitRemoval::InvSqrtE); } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1,2,3 as described above.) * @param dir The direction. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel, Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==4); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==4); } /** * Default constructor */ RSSpinorWaveFunction() : WaveFunctionBase(), _wf() {} //@} /** * Access to the wavefunction and its components. */ //@{ /** * subscript operator for the wavefunction * Set components by index. */ complex operator ()(int i, int j) const { assert( i>=0 && i<=3 && j>=0 && j<=3); return _wf(i,j); } /** * return wavefunction as LorentzRSSpinor */ const LorentzRSSpinor & wave() const {return _wf;} /** * Get first spinor component for the x vector */ complex xs1() const {return _wf.xs1();} /** * Get second spinor component for the x vector */ complex xs2() const {return _wf.xs2();} /** * Get third spinor component for the x vector */ complex xs3() const {return _wf.xs3();} /** * Get fourth spinor component for the x vector */ complex xs4() const {return _wf.xs4();} /** * Get first spinor component for the y vector */ complex ys1() const {return _wf.ys1();} /** * Get second spinor component for the y vector */ complex ys2() const {return _wf.ys2();} /** * Get third spinor component for the y vector */ complex ys3() const {return _wf.ys3();} /** * Get fourth spinor component for the y vector */ complex ys4() const {return _wf.ys4();} /** * Get first spinor component for the z vector */ complex zs1() const {return _wf.zs1();} /** * Get second spinor component for the z vector */ complex zs2() const {return _wf.zs2();} /** * Get third spinor component for the z vector */ complex zs3() const {return _wf.zs3();} /** * Get fourth spinor component for the z vector */ complex zs4() const {return _wf.zs4();} /** * Get first spinor component for the t vector */ complex ts1() const {return _wf.ts1();} /** * Get second spinor component for the t vector */ complex ts2() const {return _wf.ts2();} /** * Get third spinor component for the t vector */ complex ts3() const {return _wf.ts3();} /** * Get fourth spinor component for the t vector */ complex ts4() const {return _wf.ts4();} //@} /** * reset functions */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1,2,3 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1,2,3 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * storage of the Lorentz RSSpinor */ LorentzRSSpinor _wf; /// Return wavefunction as LorentzRSSpinor LorentzRSSpinor dimensionedWf() const { LorentzRSSpinor temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) for (unsigned int j=0; j<4; ++j) temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; return temp; } }; } } #endif /* ThePEG_RSSpinorWaveFunction_H */ diff --git a/Helicity/WaveFunction/SpinorBarWaveFunction.h b/Helicity/WaveFunction/SpinorBarWaveFunction.h --- a/Helicity/WaveFunction/SpinorBarWaveFunction.h +++ b/Helicity/WaveFunction/SpinorBarWaveFunction.h @@ -1,301 +1,301 @@ // -*- C++ -*- // // SpinorBarWaveFunction.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_SpinorBarWaveFunction_H #define ThePEG_SpinorBarWaveFunction_H // // This is the declaration of the SpinorBarWaveFunction class. #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** * Forward declaration of the SpinorWaveFunction class */ class SpinorWaveFunction; /** \ingroup Helicity * \author Peter Richardson * * The SpinorBarWaveFunction class is designed to store the wavefunction * of a barred spinor in a form suitable for use in helicity amplitude * calculations of the matrix element using a similar philosophy to the * FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzSpinorBar class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and ParticleData for the fermion. * * This class also contains the code which does the actually calculation * of the barred spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$\bar{v}\f$ spinor. * - outgoing calculates a \f$\bar{u}\f$ spinor. * * N.B. In our convention 0 is the \f$-\frac12\f$ helicity state and * 1 is the \f$+\frac12\f$ helicity state * * @see WaveFunctionBase * @see LorentzSpinorBar * @see HelicityDefinitions */ class SpinorBarWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param s1 The first component * @param s2 The second component * @param s3 The third component * @param s4 The fourth component */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex s1,complex s2, complex s3,complex s4) : WaveFunctionBase(p,part), _wf(s1,s2,s3,s4) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, const LorentzSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==2); } SpinorBarWaveFunction(const tPPtr & p, const LorentzSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==2); for (unsigned int i=0; i<4; ++i) - _wf[i]=wave[i]*UnitRemoval::InvSqrtE; + _wf[i]=Complex(wave[i]*UnitRemoval::InvSqrtE); } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1 as described above.) * @param dir The direction. */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel,Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==2); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==2); } /** * Default constructor. */ SpinorBarWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ SpinorBarWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Subscript operator for the wavefunction. */ complex operator ()(int i) const { assert(i>=0 &&i<=3); return _wf(i); } /** * Return wavefunction as LorentzSpinorBar. */ const LorentzSpinorBar & wave() const {return _wf;} /// Return wavefunction as LorentzSpinorBar LorentzSpinorBar dimensionedWave() const { return dimensionedWf(); } /** * Get the first spin component component. */ complex s1() const {return _wf.s1();} /** * Get the second spin component component. */ complex s2() const {return _wf.s2();} /** * Get the third spin component component. */ complex s3() const {return _wf.s3();} /** * Get the fourth spin component component. */ complex s4() const {return _wf.s4();} /** * Take the conjugate of the spinor \f$u_c=C\bar{u}^T\f$. This operation * transforms u-spinors to v-spinors and vice-versa and is required when * dealing with majorana particles. */ void conjugate(); /** * Return the barred spinor */ SpinorWaveFunction bar(); /** * Reset functions. */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1 as described above.) */ void calculateWaveFunction(unsigned int ihel); public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Storage of the Lorentz SpinorBar wavefunction. */ LorentzSpinorBar _wf; /// Return wavefunction as LorentzSpinorBar LorentzSpinorBar dimensionedWf() const { LorentzSpinorBar temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) temp(i) = _wf(i)*UnitRemoval::SqrtE; return temp; } }; } } #endif diff --git a/Helicity/WaveFunction/SpinorWaveFunction.h b/Helicity/WaveFunction/SpinorWaveFunction.h --- a/Helicity/WaveFunction/SpinorWaveFunction.h +++ b/Helicity/WaveFunction/SpinorWaveFunction.h @@ -1,308 +1,308 @@ // -*- C++ -*- // // SpinorWaveFunction.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_SpinorWaveFunction_H #define ThePEG_SpinorWaveFunction_H // // This is the declaration of the SpinorWaveFunction class. // #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** * Forward declaration of the SpinorBarWaveFunction class */ class SpinorBarWaveFunction; /** \ingroup Helicity * \author Peter Richardson * * The SpinorWaveFunction class is designed to store the wavefunction * of a spinor in a form suitable for use in helicity amplitude calculations * of the matrix element using a similar philosophy to the FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzSpinor class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and getParticleData for the fermion. * * This class also contains the code which does the actually calculation * of the spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$u\f$ spinor. * - outgoing calculates a \f$v\f$ spinor. * * N.B. In our convention 0 is the \f$-\frac12\f$ helicity state and * 1 is the \f$+\frac12\f$ helicity state * * @see WaveFunctionBase * @see LorentzSpinor * @see HelicityDefinitions */ class SpinorWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param s1 The first component * @param s2 The second component * @param s3 The third component * @param s4 The fourth component */ SpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part,complex s1, complex s2,complex s3,complex s4) : WaveFunctionBase(p,part), _wf(s1,s2,s3,s4) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part, const LorentzSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The particle * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorWaveFunction(const tPPtr & p, const LorentzSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(), p->dataPtr(), dir), _wf(wave.Type()) { assert(iSpin()==2); for (unsigned int i=0; i<4; ++i) - _wf[i]=wave[i]*UnitRemoval::InvSqrtE; + _wf[i]=Complex(wave[i]*UnitRemoval::InvSqrtE); } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1 as described above.) * @param dir The direction. */ SpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel, Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==2); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ SpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==2); } /** * Default constructor. */ SpinorWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ SpinorWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Subscript operator for the wavefunction. */ complex operator ()(int i) const { assert(i>=0 &&i<=3); return _wf(i); } /** * Return wavefunction as LorentzSpinor. */ const LorentzSpinor & wave() const {return _wf;} /** * Return wavefunction as LorentzSpinor. */ LorentzSpinor dimensionedWave() const {return dimensionedWf();} /** * Get the first spin component component. */ complex s1() const {return _wf.s1();} /** * Get the second spin component component. */ complex s2() const {return _wf.s2();} /** * Get the third spin component component. */ complex s3() const {return _wf.s3();} /** * Get the fourth spin component component. */ complex s4() const {return _wf.s4();} //@} /** * Take the conjugate of the spinor \f$u_c=C\bar{u}^T\f$. This operation * transforms u-spinors to v-spinors and vice-versa and is required when * dealing with majorana particles. */ void conjugate() { _wf=_wf.conjugate(); } /** * Return the barred spinor */ SpinorBarWaveFunction bar(); /** * Reset functions. */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * Storage of the Lorentz Spinor. */ LorentzSpinor _wf; /// Return wavefunction as LorentzSpinor LorentzSpinor dimensionedWf() const { LorentzSpinor temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) temp(i) = _wf(i)*UnitRemoval::SqrtE; return temp; } }; } } #endif diff --git a/Helicity/WaveFunction/TensorWaveFunction.cc b/Helicity/WaveFunction/TensorWaveFunction.cc --- a/Helicity/WaveFunction/TensorWaveFunction.cc +++ b/Helicity/WaveFunction/TensorWaveFunction.cc @@ -1,338 +1,338 @@ // -*- C++ -*- // // TensorWaveFunction.cc is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2003-2017 Peter Richardson, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // // // This is the implementation of the non-inlined, non-templated member // functions of the TensorWaveFunction class. // // Author: Peter Richardson // #include "TensorWaveFunction.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the actual wavefunction void TensorWaveFunction::calculateWaveFunction(unsigned int ihel, TensorPhase tphase) { int jhel=ihel-2; assert(direction()!=intermediate); // check for a valid helicty combination assert((jhel<=2 && jhel>=-2 && mass() >ZERO) || ((jhel==2 || jhel==-2) && mass()==ZERO)); // extract the momentum components double fact = direction()==outgoing ? -1. : 1; Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass(); // calculate some kinematic quantites; Energy2 pt2 = sqr(ppx)+sqr(ppy); Energy pabs = sqrt(pt2+sqr(ppz)); Energy pt = sqrt(pt2); // polarization vectors complex epsp[4],epsm[4],eps0[4]; // + helicity vector if needed if(jhel>=0) { // calculate the overall phase complexphase; if(tphase==tensor_phase) { phase = pt==ZERO ? 1. : complex(ppx/pt,-fact*ppy/pt); } else phase = 1.; phase = phase*sqrt(0.5); // first the no pt case if(pt==ZERO) { double sgnz = ppz(0,-fact); epsp[2]=0.; epsp[3]=0.; } else { InvEnergy opabs=1./pabs; InvEnergy opt =1./pt; epsp[0]=phase*complex(-ppz*ppx*opabs*opt, fact*ppy*opt); epsp[1]=phase*complex(-ppz*ppy*opabs*opt, -fact*ppx*opt); - epsp[2]=pt*opabs*phase; + epsp[2]=Complex(pt*opabs*phase); epsp[3]=0.; } } // - helicity vector if needed if(jhel<=0) { // calculate the overall phase complex phase; if(tphase==tensor_phase) { phase = pt==ZERO ? 1. : complex(ppx/pt,fact*ppy/pt); } else phase = 1.; phase = phase*sqrt(0.5); // first the no pt case if(pt==ZERO) { double sgnz; if(ppz(0,-fact); epsm[2]=0.; epsm[3]=0.; } else { InvEnergy opabs=1./pabs; InvEnergy opt =1./pt; epsm[0]=phase*complex(ppz*ppx*opabs*opt, fact*ppy*opt); epsm[1]=phase*complex(ppz*ppy*opabs*opt, -fact*ppx*opt); - epsm[2]=-pt*opabs*phase; + epsm[2]=-Complex(pt*opabs*phase); epsm[3]=0.; } } // 0 helicity vector if needed if(jhel<=1 && jhel>=-1) { if(pabs==ZERO) { eps0[0] = 0.; eps0[1] = 0.; eps0[2] = 1.; eps0[3] = 0.; } else { InvEnergy empabs=pee/pmm/pabs; eps0[0] = empabs*ppx; eps0[1] = empabs*ppy; eps0[2] = empabs*ppz; eps0[3] = pabs/pmm; } } // put the polarization vectors together to get the wavefunction double ort; switch (jhel) { case 2: for(int ix=0;ix<4;++ix) for(int iy=0;iy<4;++iy) _wf(ix,iy)=epsp[ix]*epsp[iy]; break; case 1: ort = sqrt(0.5); for(int ix=0;ix<4;++ix) for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsp[ix]*eps0[iy]+ eps0[ix]*epsp[iy]); break; case 0: ort = 1./sqrt(6.); for(int ix=0;ix<4;++ix) for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsp[ix]*epsm[iy] + epsm[ix]*epsp[iy] +2.*eps0[ix]*eps0[iy]); break; case -1: ort = 1./sqrt(2.); for(int ix=0;ix<4;++ix) for(int iy=0;iy<4;++iy) _wf(ix,iy)=ort*( epsm[ix]*eps0[iy]+ eps0[ix]*epsm[iy]); break; case -2: for(int ix=0;ix<4;++ix) for(int iy=0;iy<4;++iy) _wf(ix,iy)=epsm[ix]*epsm[iy]; break; default: ThePEG::Helicity::HelicityConsistencyError() << "Invalid Helicity = " << jhel << " requested for Tensor" << Exception::abortnow; break; } } void TensorWaveFunction:: calculateWaveFunctions(vector > & waves, tPPtr particle,Direction dir, bool massless, TensorPhase phase) { tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(5); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<5;++ix) waves[ix]=inspin->getProductionBasisState(ix); } else { inspin->decay(); for(unsigned int ix=0;ix<5;++ix) waves[ix]=inspin->getDecayBasisState(ix); } } else { assert(!particle->spinInfo()); TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<5;++ix) { if(massless&&ix>0&&ix<5) { waves[ix] = LorentzTensor(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } } } void TensorWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle, Direction dir, bool massless, TensorPhase phase) { tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(5); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<5;++ix) waves[ix]=TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); } else { inspin->decay(); for(unsigned int ix=0;ix<5;++ix) waves[ix]=TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); } } else { assert(!particle->spinInfo()); TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<5;++ix) { if(massless&&ix>0&&ix<5) { waves[ix] = TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } } } void TensorWaveFunction:: calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction dir,bool massless, TensorPhase phase) { tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(5); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<5;++ix) waves[ix]=inspin->getProductionBasisState(ix); rho = RhoDMatrix(PDT::Spin2); } else { inspin->decay(); for(unsigned int ix=0;ix<5;++ix) waves[ix]=inspin->getDecayBasisState(ix); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<5;++ix) { if(massless&&ix>0&&ix<5) { waves[ix] = LorentzTensor(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } rho = RhoDMatrix(PDT::Spin2); } } void TensorWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle, Direction dir, bool massless, TensorPhase phase) { tTensorSpinPtr inspin = !particle->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(5); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<5;++ix) waves[ix]=TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); rho = RhoDMatrix(PDT::Spin2); } else { inspin->decay(); for(unsigned int ix=0;ix<5;++ix) waves[ix]=TensorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); TensorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<5;++ix) { if(massless&&ix>0&&ix<5) { waves[ix] = TensorWaveFunction(particle->momentum(),particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } rho = RhoDMatrix(PDT::Spin2); } } void TensorWaveFunction:: constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time,bool ) { assert(waves.size()==5); tTensorSpinPtr inspin = !part->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<5;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix]); else inspin->setDecayState(ix,waves[ix]); } else { TensorSpinPtr temp = new_ptr(TensorSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<5;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix]); else temp->setDecayState(ix,waves[ix]); } } void TensorWaveFunction:: constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool ) { assert(waves.size()==5); tTensorSpinPtr inspin = !part->spinInfo() ? tTensorSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<5;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix].wave()); else inspin->setDecayState(ix,waves[ix].wave()); } else { TensorSpinPtr temp = new_ptr(TensorSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<5;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix].wave()); else temp->setDecayState(ix,waves[ix].wave()); } } diff --git a/Vectors/LorentzVector.h b/Vectors/LorentzVector.h --- a/Vectors/LorentzVector.h +++ b/Vectors/LorentzVector.h @@ -1,783 +1,799 @@ // -*- C++ -*- // // LorentzVector.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2006-2017 David Grellscheid, Leif Lonnblad // // ThePEG is licenced under version 3 of the GPL, see COPYING for details. // Please respect the MCnet academic guidelines, see GUIDELINES for details. // #ifndef ThePEG_LorentzVector_H #define ThePEG_LorentzVector_H /** * @file LorentzVector.h contains the LorentzVector class. Lorentz * vectors can be created with any unit type as template parameter. * All basic mathematical operations are supported, as well as a * subset of the CLHEP LorentzVector functionality. */ #include "LorentzVector.fh" #include "ThePEG/Utilities/Direction.h" #include "ThePEG/Utilities/UnitIO.h" #include "LorentzRotation.h" #include "ThreeVector.h" /// Debug helper function #ifdef NDEBUG #define ERROR_IF(condition,message) if (false) {} #else #define ERROR_IF(condition,message) \ if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror) #endif namespace ThePEG { template class LorentzVector; /** * A 4-component Lorentz vector. It can be created with any unit type * as template parameter. All basic mathematical operations are * supported, as well as a subset of the CLHEP LorentzVector * functionality. */ template class LorentzVector { private: /// Value squared typedef typename BinaryOpTraits::MulT Value2; public: /** @name Constructors. */ //@{ LorentzVector() : theX(), theY(), theZ(), theT() {} LorentzVector(Value x, Value y, Value z, Value t) : theX(x), theY(y), theZ(z), theT(t) {} LorentzVector(const ThreeVector & v, Value t) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {} template LorentzVector(const LorentzVector & v) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {} //@} /// Assignment operator template LorentzVector & operator=(const LorentzVector & b) { setX(b.x()); setY(b.y()); setZ(b.z()); setT(b.t()); return *this; } public: /// @name Component access methods. //@{ Value x() const { return theX; } Value y() const { return theY; } Value z() const { return theZ; } Value t() const { return theT; } Value e() const { return t(); } //@} /// @name Component set methods. //@{ void setX(Value x) { theX = x; } void setY(Value y) { theY = y; } void setZ(Value z) { theZ = z; } void setT(Value t) { theT = t; } void setE(Value e) { setT(e); } //@} public: /// Access to the 3-component part. ThreeVector vect() const { return ThreeVector(x(),y(),z()); } /// Cast to the 3-component part. operator ThreeVector() const { return vect(); } /// Set the 3-component part. void setVect(const ThreeVector & p) { theX = p.x(); theY = p.y(); theZ = p.z(); } public: /// The complex conjugate vector. LorentzVector conjugate() const { return LorentzVector(conj(x()),conj(y()),conj(z()),conj(t())); } /// Squared magnitude \f$x^\mu\,x_\mu=t^2 - \vec{x}^2\f$. Value2 m2() const { return (t()-z())*(t()+z()) - sqr(x()) - sqr(y()); } /// Squared magnitude with another vector Value2 m2(const LorentzVector & a) const { Value tt(a.t()+t()),zz(a.z()+z()); return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y()); } /// Magnitude (signed) \f$\pm\sqrt{|t^2 - \vec{x}^2|}\f$. Value m() const { Value2 tmp = m2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Transverse mass squared \f$t^2-z^2\f$. Value2 mt2() const { return (t()-z())*(t()+z()); } /// Transverse mass (signed) \f$\pm\sqrt{|t^2 - z^2|}\f$. Value mt() const { Value2 tmp = mt2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Squared transverse component of the spatial vector \f$x^2+y^2\f$. Value2 perp2() const { return sqr(x()) + sqr(y()); } /// Transverse component of the spatial vector \f$\pm\sqrt{x^2 + y^2}\f$. Value perp() const { return sqrt(perp2()); } /** * Squared transverse component of the spatial vector with respect to the * given axis. */ template Value2 perp2(const ThreeVector & p) const { return vect().perp2(p); } /** * Transverse component of the spatial vector with respect to the * given axis. */ template Value perp(const ThreeVector & p) const { return vect().perp(p); } /// Transverse energy squared. Value2 et2() const { Value2 pt2 = vect().perp2(); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z()); } /// Transverse energy (signed). Value et() const { Value2 etet = et2(); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// Transverse energy squared with respect to the given axis. Value2 et2(const ThreeVector & v) const { Value2 pt2 = vect().perp2(v); Value pv = vect().dot(v.unit()); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv); } /// Transverse energy with respect to the given axis (signed). Value et(const ThreeVector & v) const { Value2 etet = et2(v); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// @name Spherical coordinates for the spatial part. //@{ /// Radius squared. Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); } /// Radius. Value rho() const { return sqrt(rho2()); } /// Set new radius. void setRho(Value newRho) { Value oldRho = rho(); if (oldRho == Value()) return; double factor = newRho / oldRho; setX(x()*factor); setY(y()*factor); setZ(z()*factor); } /// Polar angle. double theta() const { assert(!(x() == Value() && y() == Value() && z() == Value())); return atan2(perp(),z()); } /// Cosine of the polar angle. double cosTheta() const { Value ptot = rho(); assert( ptot > Value() ); return z() / ptot; } /// Azimuthal angle. double phi() const { return atan2(y(),x()) ; } //@} /// Pseudorapidity of spatial part. double eta() const { Value m = rho(); if ( m == Value() ) return 0.0; Value pt = max(Constants::epsilon*m, perp()); double rap = log((m + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Spatial angle with another vector. double angle(const LorentzVector & w) const { return vect().angle(w.vect()); } /// Rapidity \f$\frac{1}{2}\ln\frac{t+z}{t-z} \f$ double rapidity() const { if ( z() == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Rapidity with respect to another vector double rapidity(const Axis & ref) const { double r = ref.mag2(); ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity"); Value vdotu = vect().dot(ref)/sqrt(r); if ( vdotu == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /** * Boost from reference frame into this vector's rest * frame: \f$\frac{\vec{x}}{t}\f$. */ Boost boostVector() const { if (t() == Value()) { if (rho2() == Value2()) return Boost(); else ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result"); } // result will make analytic sense but is physically meaningless ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector"); return vect() * (1./t()); } /** * Boost from reference frame into this vector's rest * frame: \f$-\frac{\vec{x}}{t}\f$. */ Boost findBoostToCM() const { return -boostVector(); } /// Returns the positive light-cone component \f$t + z\f$. Value plus() const { return t() + z(); } /// Returns the negative light-cone component \f$t - z\f$. Value minus() const { return t() - z(); } /// Are two vectors nearby, using Euclidean measure \f$t^2 + |\vec{x}|^2\f$? bool isNear(const LorentzVector & w, double epsilon) const { Value2 limit = abs(vect().dot(w.vect())); limit += 0.25 * sqr( t() + w.t() ); limit *= sqr(epsilon); Value2 delta = (vect() - w.vect()).mag2(); delta += sqr( t() - w.t() ); return (delta <= limit); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & transform(const SpinOneLorentzRotation & m) { return *this = m.operator*(*this); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & operator*=(const SpinOneLorentzRotation & m) { return transform(m); } /// Dot product with metric \f$(+,-,-,-)\f$ template typename BinaryOpTraits::MulT dot(const LorentzVector & a) const { return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() ); } public: /** * Apply boost. * * @param bx Component x of the boost. * @param by Component y of the boost. * @param bz Component z of the boost. * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(double bx, double by, double bz, double gamma=-1.) { const double b2 = bx*bx + by*by + bz*bz; if ( b2 == 0.0 ) return *this; if ( gamma < 0.0 ) { gamma = 1.0 / sqrt(1.0 - b2); } const Value bp = bx*x() + by*y() + bz*z(); const double gamma2 = (gamma - 1.0)/b2; setX(x() + gamma2*bp*bx + gamma*bx*t()); setY(y() + gamma2*bp*by + gamma*by*t()); setZ(z() + gamma2*bp*bz + gamma*bz*t()); setT(gamma*(t() + bp)); return *this; } /** * Apply boost. * * @param b Three-vector giving the boost. * * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(Boost b, double gamma=-1.) { return boost(b.x(), b.y(), b.z(),gamma); } /** * Apply rotation around the x-axis. * * @param phi Angle in radians. */ LorentzVector & rotateX (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value ty = y() * cosphi - z() * sinphi; theZ = z() * cosphi + y() * sinphi; theY = ty; return *this; } /** * Apply rotation around the y-axis. * * @param phi Angle in radians. */ LorentzVector & rotateY (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tz = z() * cosphi - x() * sinphi; theX = x() * cosphi + z() * sinphi; theZ = tz; return *this; } /** * Apply rotation around the z-axis. * * @param phi Angle in radians. */ LorentzVector & rotateZ (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tx = x() * cosphi - y() * sinphi; theY = y() * cosphi + x() * sinphi; theX = tx; return *this; } /** * Rotate the reference frame to a new z-axis. */ LorentzVector & rotateUz (const Axis & axis) { Axis ax = axis.unit(); double u1 = ax.x(); double u2 = ax.y(); double u3 = ax.z(); double up = u1*u1 + u2*u2; if (up>0) { up = sqrt(up); Value px = x(), py = y(), pz = z(); setX( (u1*u3*px - u2*py)/up + u1*pz ); setY( (u2*u3*px + u1*py)/up + u2*pz ); setZ( -up*px + u3*pz ); } else if (u3 < 0.) { setX(-x()); setZ(-z()); } return *this; } /** * Apply a rotation. * @param angle Rotation angle in radians. * @param axis Rotation axis. */ template LorentzVector & rotate(double angle, const ThreeVector & axis) { if (angle == 0.0) return *this; const U ll = axis.mag(); assert( ll > U() ); const double sa = sin(angle), ca = cos(angle); const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll; const Value xx = x(), yy = y(), zz = z(); setX((ca+(1-ca)*dx*dx) * xx +((1-ca)*dx*dy-sa*dz) * yy +((1-ca)*dx*dz+sa*dy) * zz ); setY(((1-ca)*dy*dx+sa*dz) * xx +(ca+(1-ca)*dy*dy) * yy +((1-ca)*dy*dz-sa*dx) * zz ); setZ(((1-ca)*dz*dx-sa*dy) * xx +((1-ca)*dz*dy+sa*dx) * yy +(ca+(1-ca)*dz*dz) * zz ); return *this; } public: /// @name Mathematical assignment operators. //@{ + LorentzVector & operator+=(const LorentzVector > & a) { + theX += a.x(); + theY += a.y(); + theZ += a.z(); + theT += a.t(); + return *this; + } + template LorentzVector & operator+=(const LorentzVector & a) { theX += a.x(); theY += a.y(); theZ += a.z(); theT += a.t(); return *this; } + LorentzVector & operator-=(const LorentzVector > & a) { + theX -= Complex(a.x()); + theY -= Complex(a.y()); + theZ -= Complex(a.z()); + theT -= Complex(a.t()); + return *this; + } + template LorentzVector & operator-=(const LorentzVector & a) { theX -= a.x(); theY -= a.y(); theZ -= a.z(); theT -= a.t(); return *this; } LorentzVector & operator*=(double a) { theX *= a; theY *= a; theZ *= a; theT *= a; return *this; } LorentzVector & operator/=(double a) { theX /= a; theY /= a; theZ /= a; theT /= a; return *this; } //@} private: /// @name Vector components //@{ Value theX; Value theY; Value theZ; Value theT; //@} }; /// @name Basic mathematical operations //@{ template inline LorentzVector operator/(const LorentzVector & v, Value a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } inline LorentzVector operator/(const LorentzVector & v, Complex a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } template inline LorentzVector operator-(const LorentzVector & v) { return LorentzVector(-v.x(),-v.y(),-v.z(),-v.t()); } template inline LorentzVector operator+(LorentzVector a, const LorentzVector & b) { return a += b; } template inline LorentzVector operator-(LorentzVector a, const LorentzVector & b) { return a -= b; } template inline LorentzVector operator*(const LorentzVector & a, double b) { return LorentzVector(a.x()*b, a.y()*b, a.z()*b, a.t()*b); } template inline LorentzVector operator*(double b, LorentzVector a) { return a *= b; } template inline LorentzVector::MulT> operator*(ValueB a, const LorentzVector & v) { typedef typename BinaryOpTraits::MulT ResultT; return LorentzVector(a*v.x(), a*v.y(), a*v.z(), a*v.t()); } template inline LorentzVector::MulT> operator*(const LorentzVector & v, ValueB b) { return b*v; } template inline LorentzVector::DivT> operator/(const LorentzVector & v, ValueB b) { typedef typename BinaryOpTraits::DivT ResultT; return LorentzVector(v.x()/b, v.y()/b, v.z()/b, v.t()/b); } //@} /// @name Scalar product with metric \f$(+,-,-,-)\f$ //@{ template inline typename BinaryOpTraits::MulT operator*(const LorentzVector & a, const LorentzVector & b) { return a.dot(b); } template inline typename BinaryOpTraits::MulT operator*(const LorentzVector & a, const LorentzVector & b) { return a.dot(b); } //@} /// Equality template inline bool operator==(const LorentzVector & a, const LorentzVector & b) { return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t(); } /// Stream output. Format \f$(x,y,z;t)\f$. inline ostream & operator<< (ostream & os, const LorentzVector & v) { return os << "(" << v.x() << "," << v.y() << "," << v.z() << ";" << v.t() << ")"; } /** Return the positive light-cone component. Or negative if the * current Direction<0> is reversed. */ template inline Value dirPlus(const LorentzVector & p) { return Direction<0>::pos()? p.plus(): p.minus(); } /** Return the negative light-cone component. Or positive if the * current Direction<0> is reversed. */ template inline Value dirMinus(const LorentzVector & p) { return Direction<0>::neg()? p.plus(): p.minus(); } /** Return the component along the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline Value dirZ(const LorentzVector & p) { return Direction<0>::dir()*p.z(); } /** Return the polar angle wrt. the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline double dirTheta(const LorentzVector & p) { return Direction<0>::pos()? p.theta(): Constants::pi - p.theta(); } /** Return the cosine of the polar angle wrt. the positive z-axis. Or * the negative z-axis if the current Direction<0> is reversed. */ template inline double dirCosTheta(const LorentzVector & p) { return Direction<0>::pos()? p.cosTheta(): -p.cosTheta(); } /** Get the boost vector for the LorentzVector. If the current * Direction<0> is reversed, so is the z-component. */ template inline ThreeVector dirBoostVector(const LorentzVector & p) { ThreeVector b(p.boostVector()); if ( Direction<0>::neg() ) b.setZ(-b.z()); return b; } /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Value x, Value y) { LorentzVector r(x, y, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone components. */ template inline LorentzVector lightCone(Value plus, Value minus) { // g++-3.3 has a problem with using Value() directly // gcc-bug c++/3650, fixed in 3.4 static const Value zero = Value(); LorentzVector r(zero, zero, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } } // delayed header inclusion to break inclusion loop: // LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec #include "Transverse.h" namespace ThePEG { /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Value x = Value(), Value y = Value()) { LorentzVector r(x, y, Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Output a LorentzVector with units to a stream. */ template void ounitstream(OStream & os, const LorentzVector & p, UnitT & u) { os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u) << ounit(p.e(), u); } /** Input a LorentzVector with units from a stream. */ template void iunitstream(IStream & is, LorentzVector & p, UnitT & u) { Value x, y, z, e; is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u); p = LorentzVector(x, y, z, e); } /// @name Traits for binary operations //@{ template struct BinaryOpTraits; /** * Template for multiplication by scalar */ template struct BinaryOpTraits, U> { /** The type resulting from multiplication of the template type with itself. */ typedef LorentzVector::MulT> MulT; /** The type resulting from division of one template type with another. */ typedef LorentzVector::DivT> DivT; }; /** * Template for multiplication by scalar */ template struct BinaryOpTraits > { /** The type resulting from multiplication of the template type with itself. */ typedef LorentzVector::MulT> MulT; }; //@} } #undef ERROR_IF #endif /* ThePEG_LorentzVector_H */