Page MenuHomeHEPForge
Files F8309517

No OneTemporary
Actions

View Options

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<ZERO ? -1. : 1.;
+ return LorentzPolarizationVector(-complex<double>(jhel)*phase,
+ sgnz*phase*complex<double>(0,-fact),
+ 0.,0.);
+ }
+ else {
+ InvEnergy opabs=1./pabs;
+ InvEnergy opt =1./pt;
+ return LorentzPolarizationVector(phase*complex<double>(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt),
+ phase*complex<double>(-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<ZERO ? -1. : 1.;
- _wf.setX(-complex<double>(jhel)*phase);
- _wf.setY(sgnz*phase*complex<double>(0,-fact));
- _wf.setZ(0.);
- _wf.setT(0.);
- }
- else {
- InvEnergy opabs=1./pabs;
- InvEnergy opt =1./pt;
- _wf.setX(phase*complex<double>(-jhel*ppz*ppx*opabs*opt, fact*ppy*opt));
- _wf.setY(phase*complex<double>(-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<LorentzPolarizationVector> & waves,
tPPtr particle,Direction dir,bool massless,
VectorPhase phase) {
tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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<VectorWaveFunction> & waves,
tPPtr particle, Direction dir, bool massless,
VectorPhase phase) {
tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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<LorentzPolarizationVector> & waves,
RhoDMatrix & rho,tPPtr particle,
Direction dir,bool massless,
VectorPhase phase) {
tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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<VectorWaveFunction> & waves,
RhoDMatrix & rho,
tPPtr particle, Direction dir,bool massless,
VectorPhase phase) {
tVectorSpinPtr inspin = !particle->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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<LorentzPolarizationVector> & waves,
tPPtr part,Direction dir, bool time,bool ) {
assert(waves.size()==3);
tVectorSpinPtr inspin = !part->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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<VectorWaveFunction> & waves,
tPPtr part,Direction dir, bool time,bool ) {
assert(waves.size()==3);
tVectorSpinPtr inspin = !part->spinInfo() ? tVectorSpinPtr() :
dynamic_ptr_cast<tVectorSpinPtr>(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());
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 3:36 PM (1 d, 8 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023311
Default Alt Text
(13 KB)

Event Timeline