diff --git a/Helicity/HelicityFunctions.h b/Helicity/HelicityFunctions.h new file mode 100644 --- /dev/null +++ b/Helicity/HelicityFunctions.h @@ -0,0 +1,82 @@ +// -*- 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" + +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 + if(ihel==10) return p*UnitRemoval::InvE; + // check a valid helicity combination + assert(ihel==0 || ihel == 2 || (ihel==1 && 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); + // 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)); + } + } +} + +} +} +} + +#endif /* ThePEG_HelicityFunctions_H */ diff --git a/Helicity/Makefile.am b/Helicity/Makefile.am --- a/Helicity/Makefile.am +++ b/Helicity/Makefile.am @@ -1,33 +1,33 @@ SUBDIRS = WaveFunction Vertex mySOURCES = LorentzSpinor.cc \ LorentzSpinorBar.cc LorentzTensor.cc \ FermionSpinInfo.cc VectorSpinInfo.cc ScalarSpinInfo.cc TensorSpinInfo.cc \ LorentzRSSpinor.cc LorentzRSSpinorBar.cc RSFermionSpinInfo.cc DOCFILES = LorentzPolarizationVector.h LorentzSpinor.h \ LorentzSpinorBar.h LorentzTensor.h \ FermionSpinInfo.h VectorSpinInfo.h ScalarSpinInfo.h TensorSpinInfo.h \ LorentzRSSpinor.h LorentzRSSpinorBar.h RSFermionSpinInfo.h \ - epsilon.h + epsilon.h HelicityFunctions.h INCLUDEFILES = $(DOCFILES) \ LorentzSpinor.fh LorentzSpinor.tcc LorentzSpinorBar.tcc \ LorentzSpinorBar.fh \ FermionSpinInfo.fh \ LorentzTensor.fh \ VectorSpinInfo.fh \ ScalarSpinInfo.fh \ TensorSpinInfo.fh LorentzTensor.tcc \ HelicityDefinitions.h \ LorentzRSSpinor.fh LorentzRSSpinorBar.fh RSFermionSpinInfo.fh \ LorentzRSSpinor.tcc LorentzRSSpinorBar.tcc noinst_LTLIBRARIES = libThePEGHelicity.la libThePEGHelicity_la_SOURCES = $(mySOURCES) $(INCLUDEFILES) libThePEGHelicity_la_LIBADD = Vertex/libThePEGVertex.la \ WaveFunction/libThePEGWaveFunction.la include $(top_srcdir)/Config/Makefile.aminclude diff --git a/Helicity/WaveFunction/VectorWaveFunction.cc b/Helicity/WaveFunction/VectorWaveFunction.cc --- a/Helicity/WaveFunction/VectorWaveFunction.cc +++ b/Helicity/WaveFunction/VectorWaveFunction.cc @@ -1,306 +1,232 @@ // -*- C++ -*- // // VectorWaveFunction.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 VectorWaveFunction class. // // Author: Peter Richardson // #include "VectorWaveFunction.h" +#include "ThePEG/Helicity/HelicityFunctions.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the Wavefunction void VectorWaveFunction::calculateWaveFunction(unsigned int ihel,VectorPhase vphase) { - Direction dir=direction(); - if(dir==intermediate) - throw ThePEG::Helicity::HelicityConsistencyError() - << "In VectorWaveFunction::calcluateWaveFunction " - << "particle must be incoming or outgoing not intermediate" - << Exception::abortnow; - // check a valid helicity combination - if(ihel==0 || ihel==2||(ihel==1&&mass()>ZERO)) { - int jhel=ihel-1; - // extract the momentum components - double fact=-1.; if(dir==incoming){fact=1.;} - Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=mass(); - // calculate some kinematic quantites; - Energy2 pt2 = ppx*ppx+ppy*ppy; - Energy pabs = sqrt(pt2+ppz*ppz); - Energy pt = sqrt(pt2); - // overall phase of the vector - Complex phase; - 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); - } - else phase = 1.; - if(ihel!=1) phase = phase/sqrt(2.); - // first the +/-1 helicity states - if(ihel!=1) { - // first the no pt case - if(pt==ZERO) { - double sgnz; - sgnz = ppz(jhel)*phase); - _wf.setY(sgnz*phase*complex(0,-fact)); - _wf.setZ(0.); - _wf.setT(0.); - } - else { - InvEnergy opabs=1./pabs; - InvEnergy opt =1./pt; - _wf.setX(phase*complex(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt)); - _wf.setY(phase*complex(-jhel*ppz*ppy*opabs*opt,-fact*ppx*opt)); - _wf.setZ(double(jhel*pt*opabs)*phase); - _wf.setT(0.); - } - } - // 0 component for massive vectors - else { - if(pabs==ZERO) { - _wf.setX(0.); - _wf.setY(0.); - _wf.setZ(1.); - _wf.setT(0.); - } - else { - InvEnergy empabs=pee/pmm/pabs; - _wf.setX(double(empabs*ppx)); - _wf.setY(double(empabs*ppy)); - _wf.setZ(double(empabs*ppz)); - _wf.setT(double(pabs/pmm)); - } - } - } - // special return the momentum as a check of gauge invariance - else if(ihel==10) { - _wf.setX(double(px()/MeV)); - _wf.setY(double(py()/MeV)); - _wf.setZ(double(pz()/MeV)); - _wf.setT(double(e()/MeV)); - } - // issue warning and return zero - else { - ThePEG::Helicity::HelicityConsistencyError() - << "Invalid Helicity = " << ihel << " requested for Vector " - << particle()->PDGName() << Exception::abortnow; - _wf.setX(0.);_wf.setY(0.);_wf.setZ(0.);_wf.setT(0.); - } + _wf = HelicityFunctions::polarizationVector(momentum(),ihel,direction(),vphase); } void VectorWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle,Direction dir,bool massless, VectorPhase phase) { tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(3); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<3;++ix) waves[ix]=inspin->getProductionBasisState(ix); } else { inspin->decay(); for(unsigned int ix=0;ix<3;++ix) waves[ix]=inspin->getDecayBasisState(ix); } } else { assert(!particle->spinInfo()); VectorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<3;++ix) { if(massless&&ix==1) { waves[ix] = LorentzPolarizationVector(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } } } void VectorWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle, Direction dir, bool massless, VectorPhase phase) { tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(3); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<3;++ix) waves[ix]=VectorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); } else { inspin->decay(); for(unsigned int ix=0;ix<3;++ix) waves[ix]=VectorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); } } else { assert(!particle->spinInfo()); VectorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<3;++ix) { if(massless&&ix==1) { waves[ix] = VectorWaveFunction(particle->momentum(), particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } } } void VectorWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho,tPPtr particle, Direction dir,bool massless, VectorPhase phase) { tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(3); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<3;++ix) waves[ix]=inspin->getProductionBasisState(ix); rho = RhoDMatrix(PDT::Spin1); } else { inspin->decay(); for(unsigned int ix=0;ix<3;++ix) waves[ix]=inspin->getDecayBasisState(ix); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); VectorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<3;++ix) { if(massless&&ix==1) { waves[ix] = LorentzPolarizationVector(); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave.wave(); } } rho = RhoDMatrix(PDT::Spin1); } } void VectorWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle, Direction dir,bool massless, VectorPhase phase) { tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(3); if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<3;++ix) waves[ix]=VectorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getProductionBasisState(ix),dir); rho = RhoDMatrix(PDT::Spin1); } else { inspin->decay(); for(unsigned int ix=0;ix<3;++ix) waves[ix]=VectorWaveFunction(particle->momentum(), particle->dataPtr(), inspin->getDecayBasisState(ix),dir); rho = inspin->rhoMatrix(); } } else { assert(!particle->spinInfo()); VectorWaveFunction wave(particle->momentum(),particle->dataPtr(),0, dir,phase); for(unsigned int ix=0;ix<3;++ix) { if(massless&&ix==1) { waves[ix] = VectorWaveFunction(particle->momentum(), particle->dataPtr(),dir); } else { if(ix!=0) wave.reset(ix); waves[ix] = wave; } } rho = RhoDMatrix(PDT::Spin1); } } void VectorWaveFunction:: constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool ) { assert(waves.size()==3); tVectorSpinPtr inspin = !part->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<3;++ix) if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix]); else inspin->setDecayState(ix,waves[ix]); } else { VectorSpinPtr temp = new_ptr(VectorSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<3;++ix) if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix]); else temp->setDecayState(ix,waves[ix]); } } void VectorWaveFunction:: constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool ) { assert(waves.size()==3); tVectorSpinPtr inspin = !part->spinInfo() ? tVectorSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<3;++ix) if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix].wave()); else inspin->setDecayState(ix,waves[ix].wave()); } else { VectorSpinPtr temp = new_ptr(VectorSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<3;++ix) if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix].wave()); else temp->setDecayState(ix,waves[ix].wave()); } }