diff --git a/Helicity/HelicityFunctions.h b/Helicity/HelicityFunctions.h --- a/Helicity/HelicityFunctions.h +++ b/Helicity/HelicityFunctions.h @@ -1,280 +1,285 @@ // -*- C++ -*- // // HelicityFunctions.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_HelicityFunctions_H #define ThePEG_HelicityFunctions_H // // This is the declaration of the HelicityFunctions header for common functions // used in helicity calculations to avoid duplication of code #include "ThePEG/Vectors/LorentzVector.h" #include "LorentzSpinor.h" #include "LorentzSpinorBar.h" namespace ThePEG { namespace Helicity { namespace HelicityFunctions { inline LorentzPolarizationVector polarizationVector(const Lorentz5Momentum & p, unsigned int ihel, Direction dir, VectorPhase vphase=default_vector_phase) { // check the direction assert(dir!=intermediate); - // special helicity combination for guge invariance tests + // special helicity combination for gauge invariance tests if(ihel==10) return p*UnitRemoval::InvE; // check a valid helicity combination - assert(ihel==0 || ihel == 2 || (ihel==1 && p.mass()>ZERO)); + assert(ihel==0 || ihel == 2 || ((ihel==1 || ihel==3) && p.mass()>ZERO )); // convert the helicitty from 0,1,2 to -1,0,1 int jhel=ihel-1; // extract the momentum components double fact = dir==outgoing ? -1 : 1; Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass(); // calculate some kinematic quantites Energy2 pt2 = ppx*ppx+ppy*ppy; Energy pabs = sqrt(pt2+ppz*ppz); Energy pt = sqrt(pt2); + // zero subtracted + if(ihel==3) { + InvEnergy pre = pmm/pabs/(pee+pabs); + return LorentzPolarizationVector(double(pre*ppx),double(pre*ppy),double(pre*ppz),-double(pre*pabs)); + } // overall phase of the vector Complex phase(1.); if(vphase==vector_phase) { if(pt==ZERO || ihel==1) phase = 1.; else if(ihel==0) phase = Complex(ppx/pt,-fact*ppy/pt); else phase = Complex(ppx/pt, fact*ppy/pt); } if(ihel!=1) phase = phase/sqrt(2.); // first the +/-1 helicity states if(ihel!=1) { // first the zero pt case if(pt==ZERO) { double sgnz = ppz(jhel)*phase, sgnz*phase*complex(0,-fact), 0.,0.); } else { InvEnergy opabs=1./pabs; InvEnergy opt =1./pt; return LorentzPolarizationVector(phase*complex(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt), phase*complex(-jhel*ppz*ppy*opabs*opt,-fact*ppx*opt), double(jhel*pt*opabs)*phase,0.); } } // 0 component for massive vectors else { if(pabs==ZERO) { return LorentzPolarizationVector(0.,0.,1.,0.); } else { InvEnergy empabs=pee/pmm/pabs; return LorentzPolarizationVector(double(empabs*ppx),double(empabs*ppy), double(empabs*ppz),double(pabs/pmm)); } } } inline LorentzSpinor dimensionedSpinor(const Lorentz5Momentum & p, unsigned int ihel, Direction dir) { // check direction and helicity assert(dir!=intermediate); assert(ihel<=1); // extract the momentum components double fact = dir==incoming ? 1 : -1.; Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass(); // define and calculate some kinematic quantities Energy2 ptran2 = ppx*ppx+ppy*ppy; Energy pabs = sqrt(ptran2+ppz*ppz); Energy ptran = sqrt(ptran2); // first need to evalulate the 2-component helicity spinors // this is the same regardless of which definition of the spinors // we are using complex hel_wf[2]; // compute the + spinor for + helicty particles and - helicity antiparticles if((dir==incoming && ihel== 1) || (dir==outgoing && ihel==0)) { // no transverse momentum if(ptran==ZERO) { if(ppz>=ZERO) { hel_wf[0] = 1; hel_wf[1] = 0; } else { hel_wf[0] = 0; hel_wf[1] = 1; } } else { InvSqrtEnergy denominator = 1./sqrt(2.*pabs); SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); hel_wf[0] = denominator*rtppluspz; hel_wf[1] = Complex(denominator/rtppluspz*complex(ppx,ppy)); } } // compute the - spinor for - helicty particles and + helicity antiparticles else { // no transverse momentum if(ptran==ZERO) { if(ppz>=ZERO) { hel_wf[0] = 0; hel_wf[1] = 1; } // transverse momentum else { hel_wf[0] = -1; hel_wf[1] = 0; } } else { InvSqrtEnergy denominator = 1./sqrt(2.*pabs); SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); hel_wf[0] = Complex(denominator/rtppluspz*complex(-ppx,ppy)); hel_wf[1] = denominator*rtppluspz; } } SqrtEnergy upper,lower; SqrtEnergy eplusp = sqrt(max(pee+pabs,ZERO)); SqrtEnergy eminusp = ( pmm != ZERO ) ? pmm/eplusp : ZERO; // set up the coefficients for the different cases if(dir==incoming) { if(ihel==1) { upper = eminusp; lower = eplusp; } else { upper = eplusp; lower = eminusp; } } else { if(ihel==1) { upper = -eplusp; lower = eminusp; } else { upper = eminusp; lower =-eplusp; } } return LorentzSpinor(upper*hel_wf[0],upper*hel_wf[1], lower*hel_wf[0],lower*hel_wf[1], (dir==incoming) ? SpinorType::u : SpinorType::v); } inline LorentzSpinor spinor(const Lorentz5Momentum & p, unsigned int ihel, Direction dir) { LorentzSpinor temp = dimensionedSpinor(p,ihel,dir); return LorentzSpinor(temp.s1()*UnitRemoval::InvSqrtE, temp.s2()*UnitRemoval::InvSqrtE, temp.s3()*UnitRemoval::InvSqrtE, temp.s4()*UnitRemoval::InvSqrtE,temp.Type()); } inline LorentzSpinorBar dimensionedSpinorBar(const Lorentz5Momentum & p, unsigned int ihel, Direction dir) { // check direction and helicity assert(dir!=intermediate); assert(ihel<=1); // extract the momentum components double fact = dir==incoming ? 1. : -1.; Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass(); // define and calculate some kinematic quantities Energy2 ptran2 = ppx*ppx+ppy*ppy; Energy pabs = sqrt(ptran2+ppz*ppz); Energy ptran = sqrt(ptran2); // first need to evalulate the 2-component helicity spinors Complex hel_wf[2]; // compute the + spinor for + helicty particles and - helicity antiparticles if((dir==outgoing && ihel== 1) || (dir==incoming && ihel==0)) { // no transverse momentum if(ptran==ZERO) { if(ppz>=ZERO) { hel_wf[0] = 1; hel_wf[1] = 0; } else { hel_wf[0] = 0; hel_wf[1] = 1; } } else { InvSqrtEnergy denominator = 1./sqrt(2.*pabs); SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); hel_wf[0] = denominator*rtppluspz; hel_wf[1] = Complex(denominator/rtppluspz*complex(ppx,-ppy)); } } // compute the - spinor for - helicty particles and + helicity antiparticles else { // no transverse momentum if(ptran==ZERO) { if(ppz>=ZERO) { hel_wf[0] = 0; hel_wf[1] = 1; } // transverse momentum else { hel_wf[0] = -1; hel_wf[1] = 0; } } else { InvSqrtEnergy denominator = 1./sqrt(2.*pabs); SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz); hel_wf[0] = Complex(denominator/rtppluspz*complex(-ppx,-ppy)); hel_wf[1] = denominator*rtppluspz; } } SqrtEnergy upper, lower; SqrtEnergy eplusp = sqrt(max(pee+pabs,ZERO)); SqrtEnergy eminusp = ( pmm!=ZERO ) ? pmm/eplusp : ZERO; // set up the coefficients for the different cases if(dir==outgoing) { if(ihel==1) { upper = eplusp; lower = eminusp; } else { upper = eminusp; lower = eplusp; } } else { if(ihel==1) { upper = eminusp; lower = -eplusp; } else { upper =-eplusp; lower = eminusp; } } // now finally we can construct the spinors return LorentzSpinorBar(upper*hel_wf[0], upper*hel_wf[1], lower*hel_wf[0], lower*hel_wf[1], (dir==incoming) ? SpinorType::v : SpinorType::u); } inline LorentzSpinorBar spinorBar(const Lorentz5Momentum & p, unsigned int ihel, Direction dir) { LorentzSpinorBar temp = dimensionedSpinorBar(p,ihel,dir); return LorentzSpinorBar(temp.s1()*UnitRemoval::InvSqrtE, temp.s2()*UnitRemoval::InvSqrtE, temp.s3()*UnitRemoval::InvSqrtE, temp.s4()*UnitRemoval::InvSqrtE,temp.Type()); } } } } #endif /* ThePEG_HelicityFunctions_H */