diff --git a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc --- a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc +++ b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc @@ -1,330 +1,330 @@ // -*- C++ -*- // // RSSpinorBarWaveFunction.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 RSSpinorBarWaveFunction class. // #include "RSSpinorBarWaveFunction.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the Wavefunction void RSSpinorBarWaveFunction::calculateWaveFunction(unsigned int ihel) { Complex ii(0.,1.); LorentzRSSpinorBar news; if(direction()==incoming) news=LorentzRSSpinorBar(SpinorType::u); else news=LorentzRSSpinorBar(SpinorType::v); unsigned int ix,iy; assert(direction()!=intermediate); assert(ihel<=3); // only two valid helicities in massless case assert( mass()>ZERO || (ihel == 0 || ihel == 3 ) ); // extract the momentum components // compute the normal spinors to construct the RS spinor Complex hel_wf[2][2]; if(direction()==incoming) { // the + spinor hel_wf[0][0] = 0.; hel_wf[0][1] = 1.; // the - spinor hel_wf[1][0] = 1.; hel_wf[1][1] = 0.; } else { // the + spinor hel_wf[0][0] = 1.; hel_wf[0][1] = 0.; // the - spinor hel_wf[1][0] = 0.; hel_wf[1][1] = 1.; } double fact = direction()==incoming ? 1. : -1.; Energy pmm=mass(),pee=fact*e(); Energy pabs = sqrt(sqr(px())+sqr(py())+sqr(pz())); SqrtEnergy eplusp = sqrt(pee+pabs); SqrtEnergy eminusp = ( pmm == ZERO ) ? ZERO : pmm/eplusp; SqrtEnergy upper[2],lower[2]; if(direction()==incoming) { upper[0] = eminusp; lower[0] =-eplusp ; upper[1] =-eplusp ; lower[1] = eminusp; } else { upper[0] = eplusp ; lower[0] = eminusp; upper[1] = eminusp; lower[1] = eplusp ; } // now construct the spinors complex spinor[2][4]; for(ix=0;ix<2;++ix) { spinor[ix][0] = upper[ix]*hel_wf[ix][0]; spinor[ix][1] = upper[ix]*hel_wf[ix][1]; spinor[ix][2] = lower[ix]*hel_wf[ix][0]; spinor[ix][3] = lower[ix]*hel_wf[ix][1]; } // compute the polarization vectors to construct the RS spinor Complex vec[3][4]; double ort = sqrt(0.5); double r1 = ( pmm == ZERO ) ? 0. : double(pee /pmm); double r2 = ( pmm == ZERO ) ? 0. : double(pabs/pmm); if(direction()==incoming) { vec[0][0] = ort; vec[0][1] = ort*ii; vec[0][2] = 0.; vec[0][3] = 0.; vec[1][0] = 0.; vec[1][1] = 0.; vec[1][2] =-r1; vec[1][3] =-r2; vec[2][0] =-ort; vec[2][1] = ort*ii; vec[2][2] = 0.; vec[2][3] = 0.; } else { vec[0][0] =-ort; vec[0][1] = ort*ii; vec[0][2] = 0.; vec[0][3] = 0.; vec[1][0] = 0.; vec[1][1] = 0.; vec[1][2] = r1; vec[1][3] = r2; vec[2][0] = ort; vec[2][1] = ort*ii; vec[2][2] = 0.; vec[2][3] = 0.; } // now we can put the bits together to compute the RS spinor double or3(sqrt(1./3.)),tor3(sqrt(2./3.)); if(ihel==3) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy]; + news(ix,iy)=Complex(UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy]); } else if(ihel==2) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE* - (or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy]); + news(ix,iy)=Complex(UnitRemoval::InvSqrtE* + (or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy])); } else if(ihel==1) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE* - (or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy]); + news(ix,iy)=Complex(UnitRemoval::InvSqrtE* + (or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy])); } else if(ihel==0) { for(ix=0;ix<4;++ix) { for(iy=0;iy<4;++iy) { - news(ix,iy)=UnitRemoval::InvSqrtE*(vec[2][ix]*spinor[1][iy]); + news(ix,iy)=Complex(UnitRemoval::InvSqrtE*(vec[2][ix]*spinor[1][iy])); } } } // spinor is currently along the z axis, rotate so in right direction if(pabs/pmm>1e-8) { Axis axis; axis.setX(fact*momentum().x()/pabs); axis.setY(fact*momentum().y()/pabs); axis.setZ(fact*momentum().z()/pabs); double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(sinth>1e-8) { LorentzRotation rot; rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); _wf = news.transform(rot); } else if (axis.z()<0.) { LorentzRotation rot; rot.setRotateX(Constants::pi); _wf = news.transform(rot); } else { _wf = news; } } else { _wf=news; } } void RSSpinorBarWaveFunction:: calculateWaveFunctions(vector > & waves, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getProductionBasisState(ix).bar(); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getDecayBasisState(ix).bar(); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWf(); } } } void RSSpinorBarWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorBarWaveFunction(particle, inspin->getProductionBasisState(ix).bar(), dir); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).bar(), dir); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave; } } } void RSSpinorBarWaveFunction:: calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getProductionBasisState(ix).bar(); rho = RhoDMatrix(PDT::Spin3Half); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getDecayBasisState(ix).bar(); rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWf(); } rho = RhoDMatrix(PDT::Spin3Half); } } void RSSpinorBarWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorBarWaveFunction(particle, inspin->getProductionBasisState(ix).bar(), dir); rho = RhoDMatrix(PDT::Spin3Half); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorBarWaveFunction(particle, inspin->getDecayBasisState(ix).bar(), dir); rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorBarWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave; } rho = RhoDMatrix(PDT::Spin3Half); } } void RSSpinorBarWaveFunction:: constructSpinInfo(const vector > & waves, tPPtr particle,Direction dir, bool time) { assert(waves.size()==4); tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix].bar()); else inspin->setDecayState(ix,waves[ix].bar()); } else { RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time)); particle->spinInfo(temp); for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix].bar()); else temp->setDecayState(ix,waves[ix].bar()); } } void RSSpinorBarWaveFunction:: constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time) { assert(waves.size()==4); tRSFermionSpinPtr inspin = !part->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(part->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<4;++ix) if (dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf().bar()); else inspin->setDecayState(ix,waves[ix].dimensionedWf().bar()); } else { RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(part->momentum(),time)); part->spinInfo(temp); for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix].dimensionedWf().bar()); else temp->setDecayState(ix,waves[ix].dimensionedWf().bar()); } } diff --git a/Helicity/WaveFunction/RSSpinorWaveFunction.cc b/Helicity/WaveFunction/RSSpinorWaveFunction.cc --- a/Helicity/WaveFunction/RSSpinorWaveFunction.cc +++ b/Helicity/WaveFunction/RSSpinorWaveFunction.cc @@ -1,328 +1,328 @@ // -*- C++ -*- // // RSSpinorWaveFunction.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 RSSpinorWaveFunction class. // #include "RSSpinorWaveFunction.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the Wavefunction void RSSpinorWaveFunction::calculateWaveFunction(unsigned int ihel) { LorentzRSSpinor news; if(direction()==incoming) news=LorentzRSSpinor(SpinorType::u); else news=LorentzRSSpinor(SpinorType::v); unsigned int ix,iy; // check helicity and type assert(direction()!=intermediate); assert(ihel<=3); // massive // only two valid helicities in massless case assert( mass()>ZERO || ( ihel == 0 || ihel==3 ) ); // extract the momentum components // compute the normal spinors to construct the RS spinor Complex hel_wf[2][2]; if(direction()==incoming) { // the + spinor hel_wf[0][0] = 1.; hel_wf[0][1] = 0.; // the - spinor hel_wf[1][0] = 0.; hel_wf[1][1] = 1.; } else { // the + spinor hel_wf[0][0] = 0.; hel_wf[0][1] = 1.; // the - spinor hel_wf[1][0] = 1.; hel_wf[1][1] = 0.; } // prefactors double fact = direction()==incoming ? 1. : -1.; Energy pmm=mass(),pee=fact*e(); Energy pabs = sqrt(sqr(px())+sqr(py())+sqr(pz())); SqrtEnergy eplusp = sqrt(pee+pabs); SqrtEnergy eminusp = ( pmm == ZERO ) ? ZERO : pmm/eplusp; SqrtEnergy upper[2],lower[2]; if(direction()==incoming) { upper[0] = eminusp; lower[0] = eplusp ; upper[1] = eplusp ; lower[1] = eminusp; } else { upper[0] =-eplusp ; lower[0] = eminusp; upper[1] = eminusp; lower[1] =-eplusp ; } // now construct the spinors complex spinor[2][4]; for(ix=0;ix<2;++ix) { spinor[ix][0] = upper[ix]*hel_wf[ix][0]; spinor[ix][1] = upper[ix]*hel_wf[ix][1]; spinor[ix][2] = lower[ix]*hel_wf[ix][0]; spinor[ix][3] = lower[ix]*hel_wf[ix][1]; } // compute the polarization vectors to construct the RS spinor Complex vec[3][4],ii(0.,1.); double ort = sqrt(0.5); double r1 = ( pmm == ZERO ) ? 0. : double(pee /pmm); double r2 = ( pmm == ZERO ) ? 0. : double(pabs/pmm); if(direction()==incoming) { vec[0][0] =-ort; vec[0][1] =-ort*ii; vec[0][2] = 0.; vec[0][3] = 0.; vec[1][0] = 0.; vec[1][1] = 0.; vec[1][2] = r1; vec[1][3] = r2; vec[2][0] = ort; vec[2][1] =-ort*ii; vec[2][2] = 0.; vec[2][3] = 0.; } else { vec[0][0] = ort; vec[0][1] =-ort*ii; vec[0][2] = 0.; vec[0][3] = 0.; vec[1][0] = 0.; vec[1][1] = 0.; vec[1][2] =-r1; vec[1][3] =-r2; vec[2][0] =-ort; vec[2][1] =-ort*ii; vec[2][2] = 0.; vec[2][3] = 0.; } // now we can put the bits together to compute the RS spinor double or3(sqrt(1./3.)),tor3(sqrt(2./3.)); if(ihel==3) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy]; + news(ix,iy)=Complex(UnitRemoval::InvSqrtE*vec[0][ix]*spinor[0][iy]); } else if(ihel==2) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE* - (or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy]); + news(ix,iy)=Complex(UnitRemoval::InvSqrtE* + (or3*vec[0][ix]*spinor[1][iy]+tor3*vec[1][ix]*spinor[0][iy])); } else if(ihel==1) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE* - (or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy]); + news(ix,iy)=Complex(UnitRemoval::InvSqrtE* + (or3*vec[2][ix]*spinor[0][iy]+tor3*vec[1][ix]*spinor[1][iy])); } else if(ihel==0) { for(ix=0;ix<4;++ix) for(iy=0;iy<4;++iy) - news(ix,iy)=UnitRemoval::InvSqrtE*vec[2][ix]*spinor[1][iy]; + news(ix,iy)=Complex(UnitRemoval::InvSqrtE*vec[2][ix]*spinor[1][iy]); } // spinor is currently along the z axis, rotate so in right direction if(pabs/pmm>1e-8) { Axis axis; axis.setX(fact*momentum().x()/pabs); axis.setY(fact*momentum().y()/pabs); axis.setZ(fact*momentum().z()/pabs); double sinth(sqrt(sqr(axis.x())+sqr(axis.y()))); if(sinth>1e-8) { LorentzRotation rot; rot.setRotate(acos(axis.z()),Axis(-axis.y()/sinth,axis.x()/sinth,0.)); _wf= news.transform(rot); } else if (axis.z()<0.) { LorentzRotation rot; rot.setRotateX(Constants::pi); _wf= news.transform(rot); } else { _wf = news; } } else { _wf=news; } } void RSSpinorWaveFunction:: calculateWaveFunctions(vector > & waves, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getProductionBasisState(ix); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getDecayBasisState(ix); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWf(); } } } void RSSpinorWaveFunction:: calculateWaveFunctions(vector & waves, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorWaveFunction(particle, inspin->getProductionBasisState(ix),dir); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorWaveFunction(particle, inspin->getDecayBasisState(ix),dir); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave; } } } void RSSpinorWaveFunction:: calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getProductionBasisState(ix); rho = RhoDMatrix(PDT::Spin3Half); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = inspin->getDecayBasisState(ix); rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWf(); } rho = RhoDMatrix(PDT::Spin3Half); } } void RSSpinorWaveFunction:: calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction dir) { tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); waves.resize(4); // spin info object exists if(inspin) { if(dir==outgoing) { for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorWaveFunction(particle, inspin->getProductionBasisState(ix),dir); rho = RhoDMatrix(PDT::Spin3Half); } else { inspin->decay(); for(unsigned int ix=0;ix<4;++ix) waves[ix] = RSSpinorWaveFunction(particle, inspin->getDecayBasisState(ix),dir); rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); RSSpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<4;++ix) { wave.reset(ix); waves[ix] = wave; } rho = RhoDMatrix(PDT::Spin3Half); } } void RSSpinorWaveFunction:: constructSpinInfo(const vector > & waves, tPPtr particle,Direction dir,bool time) { assert(waves.size()==4); tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix]); else inspin->setDecayState(ix,waves[ix]); } else { RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time)); particle->spinInfo(temp); for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix]); else temp->setDecayState(ix,waves[ix]); } } void RSSpinorWaveFunction:: constructSpinInfo(const vector & waves, tPPtr particle,Direction dir,bool time) { assert(waves.size()==4); tRSFermionSpinPtr inspin = !particle->spinInfo() ? tRSFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) inspin->setBasisState(ix,waves[ix].dimensionedWf()); else inspin->setDecayState(ix,waves[ix].dimensionedWf()); } else { RSFermionSpinPtr temp = new_ptr(RSFermionSpinInfo(particle->momentum(),time)); particle->spinInfo(temp); for(unsigned int ix=0;ix<4;++ix) if(dir==outgoing) temp->setBasisState(ix,waves[ix].dimensionedWf()); else temp->setDecayState(ix,waves[ix].dimensionedWf()); } } 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,343 @@ // -*- 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" 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); + 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] = denominator/rtppluspz*complex(-ppx,-ppy); + 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 _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[0] = Complex(upper*hel_wf[0]*UnitRemoval::InvSqrtE); + _wf[1] = Complex(upper*hel_wf[1]*UnitRemoval::InvSqrtE); + _wf[2] = Complex(lower*hel_wf[0]*UnitRemoval::InvSqrtE); + _wf[3] = Complex(lower*hel_wf[1]*UnitRemoval::InvSqrtE); } 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()); } } } diff --git a/Helicity/WaveFunction/SpinorWaveFunction.cc b/Helicity/WaveFunction/SpinorWaveFunction.cc --- a/Helicity/WaveFunction/SpinorWaveFunction.cc +++ b/Helicity/WaveFunction/SpinorWaveFunction.cc @@ -1,340 +1,340 @@ // -*- C++ -*- // // SpinorWaveFunction.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 SpinorWaveFunction class. // // Author: Peter Richardson // #include "SpinorWaveFunction.h" #include "ThePEG/Repository/CurrentGenerator.h" #include "SpinorBarWaveFunction.h" using namespace ThePEG; using namespace Helicity; // calculate the Wavefunction void SpinorWaveFunction::calculateWaveFunction(unsigned int ihel) { // check helicity is O.K. Direction dir = direction(); if(dir==intermediate) throw ThePEG::Helicity::HelicityConsistencyError() << "In SpinorWaveFunction::calcluateWaveFunction " << "particle must be incoming or outgoing not intermediate" << Exception::abortnow; if(ihel>1) throw ThePEG::Helicity::HelicityConsistencyError() << "Invalid Helicity = " << ihel << " requested for Spinor" << Exception::abortnow; // 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(); // 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] = denominator/rtppluspz*complex(ppx,ppy); + 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] = denominator/rtppluspz*complex(-ppx,ppy); + 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; } } // now finally we can construct the spinors _wf = LorentzSpinor( (dir==incoming) ? SpinorType::u : SpinorType::v); - _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[0] = Complex(upper*hel_wf[0]*UnitRemoval::InvSqrtE); + _wf[1] = Complex(upper*hel_wf[1]*UnitRemoval::InvSqrtE); + _wf[2] = Complex(lower*hel_wf[0]*UnitRemoval::InvSqrtE); + _wf[3] = Complex(lower*hel_wf[1]*UnitRemoval::InvSqrtE); } SpinorBarWaveFunction SpinorWaveFunction::bar() { Lorentz5Momentum p = momentum(); if(direction()==outgoing) p *= -1.; tcPDPtr ptemp = particle(); if(direction()==incoming&&particle()->CC()) ptemp = particle()->CC(); return SpinorBarWaveFunction(p,ptemp,_wf.bar(),direction()); } void SpinorWaveFunction:: 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); waves[1] = inspin->getProductionBasisState(1); } 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(); waves[1] = inspin->getDecayBasisState(1).conjugate(); } else { waves[0] = inspin->getDecayBasisState(0); waves[1] = inspin->getDecayBasisState(1); } } } // do the calculation else { assert(!particle->spinInfo()); SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave.dimensionedWave(); } } } void SpinorWaveFunction:: 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] = SpinorWaveFunction(particle, inspin->getProductionBasisState(ix),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] = SpinorWaveFunction(particle, inspin->getDecayBasisState(ix).conjugate(),dir); } else { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorWaveFunction(particle, inspin->getDecayBasisState(ix),dir); } } } // do the calculation else { assert(!particle->spinInfo()); SpinorWaveFunction wave(particle->momentum(),particle->dataPtr(),dir); for(unsigned int ix=0;ix<2;++ix) { wave.reset(ix); waves[ix] = wave; } } } void SpinorWaveFunction:: 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); waves[1] = inspin->getProductionBasisState(1); 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(); waves[1] = inspin->getDecayBasisState(1).conjugate(); } else { waves[0] = inspin->getDecayBasisState(0); waves[1] = inspin->getDecayBasisState(1); } rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); SpinorWaveFunction 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 SpinorWaveFunction:: 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] = SpinorWaveFunction(particle, inspin->getProductionBasisState(ix),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] = SpinorWaveFunction(particle, inspin->getDecayBasisState(ix).conjugate(),dir); } else { for(unsigned int ix=0;ix<2;++ix) waves[ix] = SpinorWaveFunction(particle, inspin->getDecayBasisState(ix),dir); } rho = inspin->rhoMatrix(); } } // do the calculation else { assert(!particle->spinInfo()); SpinorWaveFunction 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 SpinorWaveFunction:: constructSpinInfo(const vector > & waves, tPPtr particle,Direction dir,bool time) { assert(waves.size()==2); tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix]); else inspin->setDecayState(ix,waves[ix]); } } else { FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time)); particle->spinInfo(temp); for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix]); else temp->setDecayState(ix,waves[ix]); } } } void SpinorWaveFunction:: constructSpinInfo(const vector & waves, tPPtr particle,Direction dir,bool time) { assert(waves.size()==2); tFermionSpinPtr inspin = !particle->spinInfo() ? tFermionSpinPtr() : dynamic_ptr_cast(particle->spinInfo()); if(inspin) { for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) inspin->setBasisState(ix,waves[ix].dimensionedWf()); else inspin->setDecayState(ix,waves[ix].dimensionedWf()); } } else { FermionSpinPtr temp = new_ptr(FermionSpinInfo(particle->momentum(),time)); particle->spinInfo(temp); for(unsigned int ix=0;ix<2;++ix) { if(( dir == outgoing && time) || ( dir == incoming && !time)) temp->setBasisState(ix,waves[ix].dimensionedWf()); else temp->setDecayState(ix,waves[ix].dimensionedWf()); } } } diff --git a/Vectors/LorentzVector.h b/Vectors/LorentzVector.h --- a/Vectors/LorentzVector.h +++ b/Vectors/LorentzVector.h @@ -1,799 +1,799 @@ // -*- C++ -*- // // LorentzVector.h is a part of ThePEG - Toolkit for HEP Event Generation // Copyright (C) 2006-2017 David Grellscheid, 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_LorentzVector_H #define ThePEG_LorentzVector_H /** * @file LorentzVector.h contains the LorentzVector class. Lorentz * vectors can be created with any unit type as template parameter. * All basic mathematical operations are supported, as well as a * subset of the CLHEP LorentzVector functionality. */ #include "LorentzVector.fh" #include "ThePEG/Utilities/Direction.h" #include "ThePEG/Utilities/UnitIO.h" #include "LorentzRotation.h" #include "ThreeVector.h" /// Debug helper function #ifdef NDEBUG #define ERROR_IF(condition,message) if (false) {} #else #define ERROR_IF(condition,message) \ if ( condition ) throw ThePEG::Exception( (message) , ThePEG::Exception::eventerror) #endif namespace ThePEG { template class LorentzVector; /** * A 4-component Lorentz vector. It can be created with any unit type * as template parameter. All basic mathematical operations are * supported, as well as a subset of the CLHEP LorentzVector * functionality. */ template class LorentzVector { private: /// Value squared typedef typename BinaryOpTraits::MulT Value2; public: /** @name Constructors. */ //@{ LorentzVector() : theX(), theY(), theZ(), theT() {} LorentzVector(Value x, Value y, Value z, Value t) : theX(x), theY(y), theZ(z), theT(t) {} LorentzVector(const ThreeVector & v, Value t) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(t) {} template LorentzVector(const LorentzVector & v) : theX(v.x()), theY(v.y()), theZ(v.z()), theT(v.t()) {} //@} /// Assignment operator template LorentzVector & operator=(const LorentzVector & b) { setX(b.x()); setY(b.y()); setZ(b.z()); setT(b.t()); return *this; } public: /// @name Component access methods. //@{ Value x() const { return theX; } Value y() const { return theY; } Value z() const { return theZ; } Value t() const { return theT; } Value e() const { return t(); } //@} /// @name Component set methods. //@{ void setX(Value x) { theX = x; } void setY(Value y) { theY = y; } void setZ(Value z) { theZ = z; } void setT(Value t) { theT = t; } void setE(Value e) { setT(e); } //@} public: /// Access to the 3-component part. ThreeVector vect() const { return ThreeVector(x(),y(),z()); } /// Cast to the 3-component part. operator ThreeVector() const { return vect(); } /// Set the 3-component part. void setVect(const ThreeVector & p) { theX = p.x(); theY = p.y(); theZ = p.z(); } public: /// The complex conjugate vector. LorentzVector conjugate() const { return LorentzVector(conj(x()),conj(y()),conj(z()),conj(t())); } /// Squared magnitude \f$x^\mu\,x_\mu=t^2 - \vec{x}^2\f$. Value2 m2() const { return (t()-z())*(t()+z()) - sqr(x()) - sqr(y()); } /// Squared magnitude with another vector Value2 m2(const LorentzVector & a) const { Value tt(a.t()+t()),zz(a.z()+z()); return (tt-zz)*(tt+zz)-sqr(a.x()+x())-sqr(a.y()+y()); } /// Magnitude (signed) \f$\pm\sqrt{|t^2 - \vec{x}^2|}\f$. Value m() const { Value2 tmp = m2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Transverse mass squared \f$t^2-z^2\f$. Value2 mt2() const { return (t()-z())*(t()+z()); } /// Transverse mass (signed) \f$\pm\sqrt{|t^2 - z^2|}\f$. Value mt() const { Value2 tmp = mt2(); return tmp < Value2() ? -Value(sqrt(-tmp)) : Value(sqrt(tmp)); } /// Squared transverse component of the spatial vector \f$x^2+y^2\f$. Value2 perp2() const { return sqr(x()) + sqr(y()); } /// Transverse component of the spatial vector \f$\pm\sqrt{x^2 + y^2}\f$. Value perp() const { return sqrt(perp2()); } /** * Squared transverse component of the spatial vector with respect to the * given axis. */ template Value2 perp2(const ThreeVector & p) const { return vect().perp2(p); } /** * Transverse component of the spatial vector with respect to the * given axis. */ template Value perp(const ThreeVector & p) const { return vect().perp(p); } /// Transverse energy squared. Value2 et2() const { Value2 pt2 = vect().perp2(); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+z()*z()); } /// Transverse energy (signed). Value et() const { Value2 etet = et2(); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// Transverse energy squared with respect to the given axis. Value2 et2(const ThreeVector & v) const { Value2 pt2 = vect().perp2(v); Value pv = vect().dot(v.unit()); return pt2 == Value2() ? Value2() : e()*e() * pt2/(pt2+pv*pv); } /// Transverse energy with respect to the given axis (signed). Value et(const ThreeVector & v) const { Value2 etet = et2(v); return e() < Value() ? -sqrt(etet) : sqrt(etet); } /// @name Spherical coordinates for the spatial part. //@{ /// Radius squared. Value2 rho2() const { return sqr(x()) + sqr(y()) + sqr(z()); } /// Radius. Value rho() const { return sqrt(rho2()); } /// Set new radius. void setRho(Value newRho) { Value oldRho = rho(); if (oldRho == Value()) return; double factor = newRho / oldRho; setX(x()*factor); setY(y()*factor); setZ(z()*factor); } /// Polar angle. double theta() const { assert(!(x() == Value() && y() == Value() && z() == Value())); return atan2(perp(),z()); } /// Cosine of the polar angle. double cosTheta() const { Value ptot = rho(); assert( ptot > Value() ); return z() / ptot; } /// Azimuthal angle. double phi() const { return atan2(y(),x()) ; } //@} /// Pseudorapidity of spatial part. double eta() const { Value m = rho(); if ( m == Value() ) return 0.0; Value pt = max(Constants::epsilon*m, perp()); double rap = log((m + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Spatial angle with another vector. double angle(const LorentzVector & w) const { return vect().angle(w.vect()); } /// Rapidity \f$\frac{1}{2}\ln\frac{t+z}{t-z} \f$ double rapidity() const { if ( z() == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2() + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /// Rapidity with respect to another vector double rapidity(const Axis & ref) const { double r = ref.mag2(); ERROR_IF(r == 0,"A zero vector used as reference to LorentzVector rapidity"); Value vdotu = vect().dot(ref)/sqrt(r); if ( vdotu == ZERO ) return 0.0; ERROR_IF(t() <= ZERO, "Tried to take rapidity of negative-energy Lorentz vector"); Value pt = sqrt(max(sqr(t()*Constants::epsilon), perp2(ref) + m2())); double rap = log((t() + abs(z()))/pt); return z() > ZERO? rap: -rap; } /** * Boost from reference frame into this vector's rest * frame: \f$\frac{\vec{x}}{t}\f$. */ Boost boostVector() const { if (t() == Value()) { if (rho2() == Value2()) return Boost(); else ERROR_IF(true,"boostVector computed for LorentzVector with t=0 -- infinite result"); } // result will make analytic sense but is physically meaningless ERROR_IF(m2() <= Value2(),"boostVector computed for a non-timelike LorentzVector"); return vect() * (1./t()); } /** * Boost from reference frame into this vector's rest * frame: \f$-\frac{\vec{x}}{t}\f$. */ Boost findBoostToCM() const { return -boostVector(); } /// Returns the positive light-cone component \f$t + z\f$. Value plus() const { return t() + z(); } /// Returns the negative light-cone component \f$t - z\f$. Value minus() const { return t() - z(); } /// Are two vectors nearby, using Euclidean measure \f$t^2 + |\vec{x}|^2\f$? bool isNear(const LorentzVector & w, double epsilon) const { Value2 limit = abs(vect().dot(w.vect())); limit += 0.25 * sqr( t() + w.t() ); limit *= sqr(epsilon); Value2 delta = (vect() - w.vect()).mag2(); delta += sqr( t() - w.t() ); return (delta <= limit); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & transform(const SpinOneLorentzRotation & m) { return *this = m.operator*(*this); } /// Rotate the vector. Resets \f$x^\mu\rightarrow\mathsf{M}^\mu_\nu x^\nu\f$. LorentzVector & operator*=(const SpinOneLorentzRotation & m) { return transform(m); } /// Dot product with metric \f$(+,-,-,-)\f$ template typename BinaryOpTraits::MulT dot(const LorentzVector & a) const { return t() * a.t() - ( x() * a.x() + y() * a.y() + z() * a.z() ); } public: /** * Apply boost. * * @param bx Component x of the boost. * @param by Component y of the boost. * @param bz Component z of the boost. * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(double bx, double by, double bz, double gamma=-1.) { const double b2 = bx*bx + by*by + bz*bz; if ( b2 == 0.0 ) return *this; if ( gamma < 0.0 ) { gamma = 1.0 / sqrt(1.0 - b2); } const Value bp = bx*x() + by*y() + bz*z(); const double gamma2 = (gamma - 1.0)/b2; setX(x() + gamma2*bp*bx + gamma*bx*t()); setY(y() + gamma2*bp*by + gamma*by*t()); setZ(z() + gamma2*bp*bz + gamma*bz*t()); setT(gamma*(t() + bp)); return *this; } /** * Apply boost. * * @param b Three-vector giving the boost. * * @param gamma Optional gamma parameter for higher numerical * accuracy. The user has to ensure consistency. If not given, it * will be calculated as \f$\gamma=1/\sqrt{1-\beta^2}\f$. * */ LorentzVector & boost(Boost b, double gamma=-1.) { return boost(b.x(), b.y(), b.z(),gamma); } /** * Apply rotation around the x-axis. * * @param phi Angle in radians. */ LorentzVector & rotateX (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value ty = y() * cosphi - z() * sinphi; theZ = z() * cosphi + y() * sinphi; theY = ty; return *this; } /** * Apply rotation around the y-axis. * * @param phi Angle in radians. */ LorentzVector & rotateY (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tz = z() * cosphi - x() * sinphi; theX = x() * cosphi + z() * sinphi; theZ = tz; return *this; } /** * Apply rotation around the z-axis. * * @param phi Angle in radians. */ LorentzVector & rotateZ (double phi) { double sinphi = sin(phi); double cosphi = cos(phi); Value tx = x() * cosphi - y() * sinphi; theY = y() * cosphi + x() * sinphi; theX = tx; return *this; } /** * Rotate the reference frame to a new z-axis. */ LorentzVector & rotateUz (const Axis & axis) { Axis ax = axis.unit(); double u1 = ax.x(); double u2 = ax.y(); double u3 = ax.z(); double up = u1*u1 + u2*u2; if (up>0) { up = sqrt(up); Value px = x(), py = y(), pz = z(); setX( (u1*u3*px - u2*py)/up + u1*pz ); setY( (u2*u3*px + u1*py)/up + u2*pz ); setZ( -up*px + u3*pz ); } else if (u3 < 0.) { setX(-x()); setZ(-z()); } return *this; } /** * Apply a rotation. * @param angle Rotation angle in radians. * @param axis Rotation axis. */ template LorentzVector & rotate(double angle, const ThreeVector & axis) { if (angle == 0.0) return *this; const U ll = axis.mag(); assert( ll > U() ); const double sa = sin(angle), ca = cos(angle); const double dx = axis.x()/ll, dy = axis.y()/ll, dz = axis.z()/ll; const Value xx = x(), yy = y(), zz = z(); setX((ca+(1-ca)*dx*dx) * xx +((1-ca)*dx*dy-sa*dz) * yy +((1-ca)*dx*dz+sa*dy) * zz ); setY(((1-ca)*dy*dx+sa*dz) * xx +(ca+(1-ca)*dy*dy) * yy +((1-ca)*dy*dz-sa*dx) * zz ); setZ(((1-ca)*dz*dx-sa*dy) * xx +((1-ca)*dz*dy+sa*dx) * yy +(ca+(1-ca)*dz*dz) * zz ); return *this; } public: /// @name Mathematical assignment operators. //@{ - LorentzVector & operator+=(const LorentzVector > & a) { + LorentzVector & operator+=(const LorentzVector > > & a) { theX += a.x(); theY += a.y(); theZ += a.z(); theT += a.t(); return *this; } template LorentzVector & operator+=(const LorentzVector & a) { theX += a.x(); theY += a.y(); theZ += a.z(); theT += a.t(); return *this; } - LorentzVector & operator-=(const LorentzVector > & a) { + LorentzVector & operator-=(const LorentzVector > > & a) { theX -= Complex(a.x()); theY -= Complex(a.y()); theZ -= Complex(a.z()); theT -= Complex(a.t()); return *this; } template LorentzVector & operator-=(const LorentzVector & a) { theX -= a.x(); theY -= a.y(); theZ -= a.z(); theT -= a.t(); return *this; } LorentzVector & operator*=(double a) { theX *= a; theY *= a; theZ *= a; theT *= a; return *this; } LorentzVector & operator/=(double a) { theX /= a; theY /= a; theZ /= a; theT /= a; return *this; } //@} private: /// @name Vector components //@{ Value theX; Value theY; Value theZ; Value theT; //@} }; /// @name Basic mathematical operations //@{ template inline LorentzVector operator/(const LorentzVector & v, Value a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } inline LorentzVector operator/(const LorentzVector & v, Complex a) { return LorentzVector(v.x()/a, v.y()/a, v.z()/a, v.t()/a); } template inline LorentzVector operator-(const LorentzVector & v) { return LorentzVector(-v.x(),-v.y(),-v.z(),-v.t()); } template inline LorentzVector operator+(LorentzVector a, const LorentzVector & b) { return a += b; } template inline LorentzVector operator-(LorentzVector a, const LorentzVector & b) { return a -= b; } template inline LorentzVector operator*(const LorentzVector & a, double b) { return LorentzVector(a.x()*b, a.y()*b, a.z()*b, a.t()*b); } template inline LorentzVector operator*(double b, LorentzVector a) { return a *= b; } template inline LorentzVector::MulT> operator*(ValueB a, const LorentzVector & v) { typedef typename BinaryOpTraits::MulT ResultT; return LorentzVector(a*v.x(), a*v.y(), a*v.z(), a*v.t()); } template inline LorentzVector::MulT> operator*(const LorentzVector & v, ValueB b) { return b*v; } template inline LorentzVector::DivT> operator/(const LorentzVector & v, ValueB b) { typedef typename BinaryOpTraits::DivT ResultT; return LorentzVector(v.x()/b, v.y()/b, v.z()/b, v.t()/b); } //@} /// @name Scalar product with metric \f$(+,-,-,-)\f$ //@{ template inline typename BinaryOpTraits::MulT operator*(const LorentzVector & a, const LorentzVector & b) { return a.dot(b); } template inline typename BinaryOpTraits::MulT operator*(const LorentzVector & a, const LorentzVector & b) { return a.dot(b); } //@} /// Equality template inline bool operator==(const LorentzVector & a, const LorentzVector & b) { return a.x() == b.x() && a.y() == b.y() && a.z() == b.z() && a.t() == b.t(); } /// Stream output. Format \f$(x,y,z;t)\f$. inline ostream & operator<< (ostream & os, const LorentzVector & v) { return os << "(" << v.x() << "," << v.y() << "," << v.z() << ";" << v.t() << ")"; } /** Return the positive light-cone component. Or negative if the * current Direction<0> is reversed. */ template inline Value dirPlus(const LorentzVector & p) { return Direction<0>::pos()? p.plus(): p.minus(); } /** Return the negative light-cone component. Or positive if the * current Direction<0> is reversed. */ template inline Value dirMinus(const LorentzVector & p) { return Direction<0>::neg()? p.plus(): p.minus(); } /** Return the component along the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline Value dirZ(const LorentzVector & p) { return Direction<0>::dir()*p.z(); } /** Return the polar angle wrt. the positive z-axis. Or the negative * z-axis if the current Direction<0> is reversed. */ template inline double dirTheta(const LorentzVector & p) { return Direction<0>::pos()? p.theta(): Constants::pi - p.theta(); } /** Return the cosine of the polar angle wrt. the positive z-axis. Or * the negative z-axis if the current Direction<0> is reversed. */ template inline double dirCosTheta(const LorentzVector & p) { return Direction<0>::pos()? p.cosTheta(): -p.cosTheta(); } /** Get the boost vector for the LorentzVector. If the current * Direction<0> is reversed, so is the z-component. */ template inline ThreeVector dirBoostVector(const LorentzVector & p) { ThreeVector b(p.boostVector()); if ( Direction<0>::neg() ) b.setZ(-b.z()); return b; } /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Value x, Value y) { LorentzVector r(x, y, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone components. */ template inline LorentzVector lightCone(Value plus, Value minus) { // g++-3.3 has a problem with using Value() directly // gcc-bug c++/3650, fixed in 3.4 static const Value zero = Value(); LorentzVector r(zero, zero, 0.5*(plus-minus), 0.5*(plus+minus)); return r; } } // delayed header inclusion to break inclusion loop: // LorentzVec -> Transverse -> Lorentz5Vec -> LorentzVec #include "Transverse.h" namespace ThePEG { /** Create a LorentzVector giving its light-cone and transverse * components. */ template inline LorentzVector lightCone(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), 0.5*(plus-minus), 0.5*(plus+minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Value x = Value(), Value y = Value()) { LorentzVector r(x, y, Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Create a LorentzVector giving its light-cone and transverse * components. If the current Direction<0> is reversed, so is the * z-component. */ template inline LorentzVector lightConeDir(Value plus, Value minus, Transverse pt) { LorentzVector r(pt.x(), pt.y(), Direction<0>::dir()*0.5*(plus - minus), 0.5*(plus + minus)); return r; } /** Output a LorentzVector with units to a stream. */ template void ounitstream(OStream & os, const LorentzVector & p, UnitT & u) { os << ounit(p.x(), u) << ounit(p.y(), u) << ounit(p.z(), u) << ounit(p.e(), u); } /** Input a LorentzVector with units from a stream. */ template void iunitstream(IStream & is, LorentzVector & p, UnitT & u) { Value x, y, z, e; is >> iunit(x, u) >> iunit(y, u) >> iunit(z, u) >> iunit(e, u); p = LorentzVector(x, y, z, e); } /// @name Traits for binary operations //@{ template struct BinaryOpTraits; /** * Template for multiplication by scalar */ template struct BinaryOpTraits, U> { /** The type resulting from multiplication of the template type with itself. */ typedef LorentzVector::MulT> MulT; /** The type resulting from division of one template type with another. */ typedef LorentzVector::DivT> DivT; }; /** * Template for multiplication by scalar */ template struct BinaryOpTraits > { /** The type resulting from multiplication of the template type with itself. */ typedef LorentzVector::MulT> MulT; }; //@} } #undef ERROR_IF #endif /* ThePEG_LorentzVector_H */