ThePEG/Helicity/HelicityFunctions.h

0001 // -*- C++ -*-
0002 //
0003 // HelicityFunctions.h is a part of ThePEG - Toolkit for HEP Event Generation
0004 // Copyright (C) 2003-2019 Peter Richardson, Leif Lonnblad
0005 //
0006 // ThePEG is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef ThePEG_HelicityFunctions_H
0010 #define ThePEG_HelicityFunctions_H
0011 //
0012 // This is the declaration of the HelicityFunctions header for common functions
0013 // used in helicity calculations to avoid duplication of code
0014
0015 #include "ThePEG/Vectors/LorentzVector.h"
0016 #include "LorentzSpinor.h"
0017 #include "LorentzSpinorBar.h"
0018
0019 namespace ThePEG {
0020 namespace Helicity {
0021 namespace HelicityFunctions {
0022
0023 inline LorentzPolarizationVector polarizationVector(const Lorentz5Momentum & p,
0024                             unsigned int ihel,
0025                             Direction dir,
0026                             VectorPhase vphase=default_vector_phase) {
0027   // check the direction
0028   assert(dir!=intermediate);
0029   // special helicity combination for gauge invariance tests
0030   if(ihel==10) return p*UnitRemoval::InvE;
0031   // check a valid helicity combination
0032   assert(ihel==0 || ihel == 2 || ((ihel==1 || ihel==3) && p.mass()>ZERO ));
0033   // convert the helicitty from 0,1,2 to -1,0,1
0034   int jhel=ihel-1;
0035   // extract the momentum components
0036   double fact = dir==outgoing ? -1 : 1;
0037   Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
0038   // calculate some kinematic quantites
0039   Energy2 pt2 = ppx*ppx+ppy*ppy;
0040   Energy pabs = sqrt(pt2+ppz*ppz);
0041   Energy pt = sqrt(pt2);
0042   // zero subtracted
0043   if(ihel==3) {
0044     InvEnergy pre = pmm/pabs/(pee+pabs);
0045     return LorentzPolarizationVector(double(pre*ppx),double(pre*ppy),double(pre*ppz),-double(pre*pabs));
0046   }
0047   // overall phase of the vector
0048   Complex phase(1.);
0049   if(vphase==vector_phase) {
0050     if(pt==ZERO || ihel==1) phase = 1.;
0051     else if(ihel==0)        phase = Complex(ppx/pt,-fact*ppy/pt);
0052     else                    phase = Complex(ppx/pt, fact*ppy/pt);
0053   }
0054   if(ihel!=1) phase = phase/sqrt(2.);
0055   // first the +/-1 helicity states
0056   if(ihel!=1) {
0057     // first the zero pt case
0058     if(pt==ZERO) {
0059       double sgnz = ppz<ZERO ? -1. : 1.;
0060       return LorentzPolarizationVector(-complex<double>(jhel)*phase,
0061                        sgnz*phase*complex<double>(0,-fact),
0062                        0.,0.);
0063     }
0064     else {
0065       InvEnergy opabs=1./pabs;
0066       InvEnergy opt  =1./pt;
0067       return LorentzPolarizationVector(phase*complex<double>(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt),
0068                        phase*complex<double>(-jhel*ppz*ppy*opabs*opt,-fact*ppx*opt),
0069                        double(jhel*pt*opabs)*phase,0.);
0070     }
0071   }
0072   // 0 component for massive vectors
0073   else {
0074     if(pabs==ZERO) {
0075       return LorentzPolarizationVector(0.,0.,1.,0.);
0076     }
0077     else {
0078       InvEnergy empabs=pee/pmm/pabs;
0079       return LorentzPolarizationVector(double(empabs*ppx),double(empabs*ppy),
0080                        double(empabs*ppz),double(pabs/pmm));
0081     }
0082   }
0083 }
0084
0085
0086 inline LorentzSpinor<SqrtEnergy> dimensionedSpinor(const Lorentz5Momentum & p,
0087                            unsigned int ihel,
0088                            Direction dir) {
0089   // check direction and helicity
0090   assert(dir!=intermediate);
0091   assert(ihel<=1);
0092   // extract the momentum components
0093   double fact = dir==incoming ? 1 : -1.;
0094   Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
0095   // define and calculate some kinematic quantities
0096   Energy2 ptran2  = ppx*ppx+ppy*ppy;
0097   Energy pabs   = sqrt(ptran2+ppz*ppz);
0098   Energy ptran  = sqrt(ptran2);
0099   // first need to evalulate the 2-component helicity spinors
0100   // this is the same regardless of which definition of the spinors
0101   // we are using
0102   complex <double> hel_wf[2];
0103   // compute the + spinor for + helicty particles and - helicity antiparticles
0104   if((dir==incoming && ihel== 1) || (dir==outgoing && ihel==0)) {
0105     // no transverse momentum
0106     if(ptran==ZERO) {
0107       if(ppz>=ZERO) {
0108     hel_wf[0] = 1;
0109     hel_wf[1] = 0;
0110       }
0111       else {
0112     hel_wf[0] = 0;
0113     hel_wf[1] = 1;
0114       }
0115     }
0116     else {
0117       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
0118       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
0119       hel_wf[0] = denominator*rtppluspz;
0120       hel_wf[1] = Complex(denominator/rtppluspz*complex<Energy>(ppx,ppy));
0121     }
0122   }
0123   // compute the - spinor for - helicty particles and + helicity antiparticles
0124   else {
0125     // no transverse momentum
0126     if(ptran==ZERO) {
0127       if(ppz>=ZERO) {
0128     hel_wf[0] = 0;
0129     hel_wf[1] = 1;
0130       }
0131       // transverse momentum
0132       else {
0133     hel_wf[0] = -1;
0134     hel_wf[1] =  0;
0135       }
0136     }
0137     else {
0138       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
0139       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
0140       hel_wf[0] = Complex(denominator/rtppluspz*complex<Energy>(-ppx,ppy));
0141       hel_wf[1] = denominator*rtppluspz;
0142     }
0143   }
0144
0145   SqrtEnergy upper,lower;
0146   SqrtEnergy eplusp  = sqrt(max(pee+pabs,ZERO));
0147   SqrtEnergy eminusp = ( pmm != ZERO ) ? pmm/eplusp : ZERO;
0148   // set up the coefficients for the different cases
0149   if(dir==incoming) {
0150     if(ihel==1) {
0151       upper = eminusp;
0152       lower = eplusp;
0153     }
0154     else {
0155       upper = eplusp;
0156       lower = eminusp;
0157     }
0158   }
0159   else {
0160     if(ihel==1) {
0161       upper = -eplusp;
0162       lower = eminusp;
0163     }
0164     else {
0165       upper = eminusp;
0166       lower =-eplusp;
0167     }
0168   }
0169   return LorentzSpinor<SqrtEnergy>(upper*hel_wf[0],upper*hel_wf[1],
0170                    lower*hel_wf[0],lower*hel_wf[1],
0171                    (dir==incoming) ? SpinorType::u : SpinorType::v);
0172 }
0173
0174 inline LorentzSpinor<double> spinor(const Lorentz5Momentum & p,
0175                     unsigned int ihel,
0176                     Direction dir) {
0177   LorentzSpinor<SqrtEnergy> temp = dimensionedSpinor(p,ihel,dir);
0178   return LorentzSpinor<double>(temp.s1()*UnitRemoval::InvSqrtE,
0179                    temp.s2()*UnitRemoval::InvSqrtE,
0180                    temp.s3()*UnitRemoval::InvSqrtE,
0181                    temp.s4()*UnitRemoval::InvSqrtE,temp.Type());
0182 }
0183
0184 inline LorentzSpinorBar<SqrtEnergy> dimensionedSpinorBar(const Lorentz5Momentum & p,
0185                              unsigned int ihel,
0186                              Direction dir) {
0187   // check direction and helicity
0188   assert(dir!=intermediate);
0189   assert(ihel<=1);
0190   // extract the momentum components
0191   double fact = dir==incoming ? 1. : -1.;
0192   Energy ppx=fact*p.x(),ppy=fact*p.y(),ppz=fact*p.z(),pee=fact*p.e(),pmm=p.mass();
0193   // define and calculate some kinematic quantities
0194   Energy2 ptran2  = ppx*ppx+ppy*ppy;
0195   Energy pabs   = sqrt(ptran2+ppz*ppz);
0196   Energy ptran  = sqrt(ptran2);
0197   // first need to evalulate the 2-component helicity spinors
0198   Complex hel_wf[2];
0199   // compute the + spinor for + helicty particles and - helicity antiparticles
0200   if((dir==outgoing && ihel== 1) || (dir==incoming && ihel==0)) {
0201     // no transverse momentum
0202     if(ptran==ZERO) {
0203       if(ppz>=ZERO) {
0204     hel_wf[0] = 1;
0205     hel_wf[1] = 0;
0206       }
0207       else {
0208     hel_wf[0] = 0;
0209     hel_wf[1] = 1;
0210       }
0211     }
0212     else {
0213       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
0214       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
0215       hel_wf[0] = denominator*rtppluspz;
0216       hel_wf[1] = Complex(denominator/rtppluspz*complex<Energy>(ppx,-ppy));
0217     }
0218   }
0219   // compute the - spinor for - helicty particles and + helicity antiparticles
0220   else {
0221     // no transverse momentum
0222     if(ptran==ZERO) {
0223       if(ppz>=ZERO) {
0224     hel_wf[0] = 0;
0225     hel_wf[1] = 1;
0226       }
0227       // transverse momentum
0228       else {
0229     hel_wf[0] = -1;
0230     hel_wf[1] =  0;
0231       }
0232     }
0233     else {
0234       InvSqrtEnergy denominator = 1./sqrt(2.*pabs);
0235       SqrtEnergy rtppluspz = (ppz>=ZERO) ? sqrt(pabs+ppz) : ptran/sqrt(pabs-ppz);
0236       hel_wf[0] = Complex(denominator/rtppluspz*complex<Energy>(-ppx,-ppy));
0237       hel_wf[1] = denominator*rtppluspz;
0238     }
0239   }
0240   SqrtEnergy upper, lower;
0241   SqrtEnergy eplusp  = sqrt(max(pee+pabs,ZERO));
0242   SqrtEnergy eminusp = ( pmm!=ZERO ) ? pmm/eplusp : ZERO;
0243   // set up the coefficients for the different cases
0244   if(dir==outgoing) {
0245     if(ihel==1) {
0246       upper = eplusp;
0247       lower = eminusp;
0248     }
0249     else {
0250       upper = eminusp;
0251       lower = eplusp;
0252     }
0253   }
0254   else {
0255     if(ihel==1) {
0256     upper = eminusp;
0257     lower = -eplusp;
0258     }
0259     else {
0260       upper =-eplusp;
0261       lower = eminusp;
0262     }
0263   }
0264   // now finally we can construct the spinors
0265   return LorentzSpinorBar<SqrtEnergy>(upper*hel_wf[0],
0266                       upper*hel_wf[1],
0267                       lower*hel_wf[0],
0268                       lower*hel_wf[1],
0269                       (dir==incoming) ? SpinorType::v : SpinorType::u);
0270 }
0271
0272 inline LorentzSpinorBar<double> spinorBar(const Lorentz5Momentum & p,
0273                       unsigned int ihel,
0274                       Direction dir) {
0275   LorentzSpinorBar<SqrtEnergy> temp = dimensionedSpinorBar(p,ihel,dir);
0276   return LorentzSpinorBar<double>(temp.s1()*UnitRemoval::InvSqrtE,
0277                   temp.s2()*UnitRemoval::InvSqrtE,
0278                   temp.s3()*UnitRemoval::InvSqrtE,
0279                   temp.s4()*UnitRemoval::InvSqrtE,temp.Type());
0280 }
0281 }
0282 }
0283 }
0284
0285 #endif /* ThePEG_HelicityFunctions_H */