diff --git a/Helicity/HelicityFunctions.h b/Helicity/HelicityFunctions.h --- a/Helicity/HelicityFunctions.h +++ b/Helicity/HelicityFunctions.h @@ -1,82 +1,172 @@ // -*- 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 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); + 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 LorentzSpinorBar spinorBar(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] = 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] = 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]*UnitRemoval::InvSqrtE, + upper*hel_wf[1]*UnitRemoval::InvSqrtE, + lower*hel_wf[0]*UnitRemoval::InvSqrtE, + lower*hel_wf[1]*UnitRemoval::InvSqrtE, + (dir==incoming) ? SpinorType::v : SpinorType::u); +} + } } } #endif /* ThePEG_HelicityFunctions_H */ diff --git a/Helicity/WaveFunction/SpinorBarWaveFunction.cc b/Helicity/WaveFunction/SpinorBarWaveFunction.cc --- a/Helicity/WaveFunction/SpinorBarWaveFunction.cc +++ b/Helicity/WaveFunction/SpinorBarWaveFunction.cc @@ -1,343 +1,254 @@ // -*- C++ -*- // // SpinorBarWaveFunction.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 SpinorBarWaveFunction class. // // Author: Peter Richardson // #include "SpinorBarWaveFunction.h" #include "SpinorWaveFunction.h" +#include "ThePEG/Helicity/HelicityFunctions.h" using namespace ThePEG; using namespace Helicity; // calculate the Wavefunction void SpinorBarWaveFunction::calculateWaveFunction(unsigned int ihel) { - Direction dir=direction(); - if(dir==intermediate) ThePEG::Helicity::HelicityConsistencyError() - << "In SpinorBarWaveFunction::calcluateWaveFunction " - << "particle must be incoming or outgoing not intermediate" - << Exception::abortnow; - // check ihelicity is O.K. - if(ihel>1) ThePEG::Helicity::HelicityConsistencyError() - << "Invalid Helicity = " << ihel << " requested for SpinorBar" - << Exception::abortnow; - // extract the momentum components - double fact = dir==incoming ? 1. : -1.; - Energy ppx=fact*px(),ppy=fact*py(),ppz=fact*pz(),pee=fact*e(),pmm=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==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] = 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] = 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 - _wf = LorentzSpinorBar((dir==incoming) ? SpinorType::v : SpinorType::u); - _wf[0] = upper*hel_wf[0]*UnitRemoval::InvSqrtE; - _wf[1] = upper*hel_wf[1]*UnitRemoval::InvSqrtE; - _wf[2] = lower*hel_wf[0]*UnitRemoval::InvSqrtE; - _wf[3] = lower*hel_wf[1]*UnitRemoval::InvSqrtE; + _wf = HelicityFunctions::spinorBar(momentum(),ihel,direction()); } void SpinorBarWaveFunction::conjugate() { _wf=_wf.conjugate(); } SpinorWaveFunction SpinorBarWaveFunction::bar() { Lorentz5Momentum p = momentum(); if(direction()==outgoing) p *= -1.; tcPDPtr ptemp = particle(); if(direction()==incoming&&particle()->CC()) ptemp = particle()->CC(); return SpinorWaveFunction(p,ptemp,_wf.bar(),direction()); } void SpinorBarWaveFunction:: calculateWaveFunctions(vector > & waves, tPPtr particle,Direction dir) { tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(2); // spin info object exists if(inspin) { if(dir==outgoing) { waves[0] = inspin->getProductionBasisState(0).bar(); waves[1] = inspin->getProductionBasisState(1).bar(); } else { inspin->decay(); if( (particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) || (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) { waves[0] = inspin->getDecayBasisState(0).conjugate().bar(); waves[1] = inspin->getDecayBasisState(1).conjugate().bar(); } else { waves[0] = inspin->getDecayBasisState(0).bar(); waves[1] = inspin->getDecayBasisState(1).bar(); } } } // do the calculation else { assert(!particle->spinInfo()); SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWave(); } } } void SpinorBarWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle,Direction dir) { tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(2); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getProductionBasisState(ix).bar(), dir); } else { inspin->decay(); if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) || (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).conjugate().bar(),dir); } else { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).bar(),dir); } } } // do the calculation else { assert(!particle->spinInfo()); SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave; } } } void SpinorBarWaveFunction:: calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(2); // spin info object exists if(inspin) { if(dir==outgoing) { waves[0] = inspin->getProductionBasisState(0).bar(); waves[1] = inspin->getProductionBasisState(1).bar(); rho = RhoDMatrix(PDT::Spin1Half); } else { inspin->decay(); if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) || (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) { waves[0] = inspin->getDecayBasisState(0).conjugate().bar(); waves[1] = inspin->getDecayBasisState(1).conjugate().bar(); } else { waves[0] = inspin->getDecayBasisState(0).bar(); waves[1] = inspin->getDecayBasisState(1).bar(); } rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWave(); } rho = RhoDMatrix(PDT::Spin1Half); } } void SpinorBarWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(2); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getProductionBasisState(ix).bar(), dir); rho = RhoDMatrix(PDT::Spin1Half); } else { inspin->decay(); if((particle->id()>0&&inspin->getDecayBasisState(0).Type()!=SpinorType::u) || (particle->id()<0&&inspin->getDecayBasisState(0).Type()!=SpinorType::v)) { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).conjugate().bar(),dir); } else { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).bar(),dir); } rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); SpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave; } rho = RhoDMatrix(PDT::Spin1Half); } } void SpinorBarWaveFunction:: constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time) { assert(waves.size()==2); tFermionSpinPtr inspin = !part->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix].bar()); else inspin->setDecayState(ix,waves[ix].bar()); } } else { FermionSpinPtr temp = new_ptr(FermionSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix].bar()); else temp->setDecayState(ix,waves[ix].bar()); } } } void SpinorBarWaveFunction:: constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time) { assert(waves.size()==2); tFermionSpinPtr inspin = !part->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<2;++ix) if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix].dimensionedWf().bar()); else inspin->setDecayState(ix,waves[ix].dimensionedWf().bar()); } else { FermionSpinPtr temp = new_ptr(FermionSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix].dimensionedWf().bar()); else temp->setDecayState(ix,waves[ix].dimensionedWf().bar()); } } }