diff --git a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc --- a/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc +++ b/Helicity/WaveFunction/RSSpinorBarWaveFunction.cc @@ -1,274 +1,282 @@ // -*- 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" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the Wavefunction void RSSpinorBarWaveFunction::calculateWaveFunction(unsigned int ihel) { // if rest frame, make sure speical case is used Lorentz5Momentum ptemp=momentum(); double pr = ptemp.vect().mag2()/sqr(momentum().e()); if(pr<1e-20) { ptemp.setX(ZERO); ptemp.setY(ZERO); ptemp.setZ(ZERO); } assert(direction()!=intermediate); assert(ihel<=3); // only two valid helicities in massless case assert( mass()>ZERO || (ihel == 0 || ihel == 3 ) ); // new direct calculation _wf = LorentzRSSpinorBar(direction()==outgoing ? SpinorType::u : SpinorType::v); // compute the spinors std::array,2> spin; if(ihel!=0) spin[0] = HelicityFunctions::spinorBar(ptemp,1,direction()); if(ihel!=3) spin[1] = HelicityFunctions::spinorBar(ptemp,0,direction()); // compute the polarization vectors to construct the RS spinor std::array eps; if(ihel>=2) eps[0] = HelicityFunctions::polarizationVector(ptemp,2,direction(),vector_phase); else eps[2] = HelicityFunctions::polarizationVector(ptemp,0,direction(),vector_phase); if(mass()!=ZERO&&ihel!=0&&ihel!=3) eps[1] = HelicityFunctions::polarizationVector(ptemp,1,direction()); // 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(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = eps[0].x()*spin[0][iy]; _wf(1,iy) = eps[0].y()*spin[0][iy]; _wf(2,iy) = eps[0].z()*spin[0][iy]; _wf(3,iy) = eps[0].t()*spin[0][iy]; } } else if(ihel==2) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = or3*eps[0].x()*spin[1][iy]+tor3*eps[1].x()*spin[0][iy]; _wf(1,iy) = or3*eps[0].y()*spin[1][iy]+tor3*eps[1].y()*spin[0][iy]; _wf(2,iy) = or3*eps[0].z()*spin[1][iy]+tor3*eps[1].z()*spin[0][iy]; _wf(3,iy) = or3*eps[0].t()*spin[1][iy]+tor3*eps[1].t()*spin[0][iy]; } } else if(ihel==1) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = or3*eps[2].x()*spin[0][iy]+tor3*eps[1].x()*spin[1][iy]; _wf(1,iy) = or3*eps[2].y()*spin[0][iy]+tor3*eps[1].y()*spin[1][iy]; _wf(2,iy) = or3*eps[2].z()*spin[0][iy]+tor3*eps[1].z()*spin[1][iy]; _wf(3,iy) = or3*eps[2].t()*spin[0][iy]+tor3*eps[1].t()*spin[1][iy]; } } else if(ihel==0) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = eps[2].x()*spin[1][iy]; _wf(1,iy) = eps[2].y()*spin[1][iy]; _wf(2,iy) = eps[2].z()*spin[1][iy]; _wf(3,iy) = eps[2].t()*spin[1][iy]; } } // this makes the phase choice the same as madgraph, useful for debugging only // Energy pt = ptemp.perp(); // double fact = direction()==incoming ? 1. : -1.; // Complex emphi(fact*ptemp.x()/pt,-fact*ptemp.y()/pt); // if(ihel==3) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) /= emphi; // } // else if(ihel==1) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) *= emphi; // } // else if(ihel==0) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) *= sqr(emphi); // } } 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; + calculateWaveFunctions(waves,particle->momentum(),particle->dataPtr(),dir); } +} + +void RSSpinorBarWaveFunction:: +calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction dir) { + waves.resize(4); + RSSpinorBarWaveFunction wave(momentum,parton,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/RSSpinorBarWaveFunction.h b/Helicity/WaveFunction/RSSpinorBarWaveFunction.h --- a/Helicity/WaveFunction/RSSpinorBarWaveFunction.h +++ b/Helicity/WaveFunction/RSSpinorBarWaveFunction.h @@ -1,371 +1,378 @@ // -*- C++ -*- // // RSSpinorBarWaveFunction.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_RSSpinorBarWaveFunction_H #define ThePEG_RSSpinorBarWaveFunction_H // // This is the declaration of the RSSpinorBarWaveFunction class. // #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** \ingroup Helicity * * The RSSpinorBarWaveFunction class is designed to * store the wavefunction of a spin-\f$\frac32\f$ particle in a form * suitable for use in helicity amplitude calculations of the matrix * element using a similar philosophy to the FORTRAN HELAS code. * * In addition to storing the barred spinor using the * LorentzRSSpinorBar class it inherits from the * WaveFunctionBase class to provide storage of the * momentum and ParticleData for the fermion. * * This class also contains the code which does the actually * calculation of the barred spinor for an external particle * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$\bar{v}\f$ spinor. * - outgoing calculates a \f$\bar{u}\f$ spinor. * * The barred spinors are calculated using a Clebsch-Gordon decomposition * in the rest-frame for a massive particle and boosted to the lab-frame. * For massless particles the calculation is performed in the lab-frame * (N.B. there are only two helicities \f$\pm\frac32\f$ in this case.) * * N.B. In our convention 0 is the \f$-\frac32\f$ helicity state, * 1 is the \f$-\frac12\f$ helicity state, * 2 is the \f$+\frac12\f$ helicity state * 3 is the \f$+\frac32\f$ helicity state and * * @see WaveFunctionBase * @see LorentzRSSpinorBar * @see HelicityDefinitions * */ class RSSpinorBarWaveFunction: public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor * @param p The momentum. * @param part The ParticleData pointer. * @param xs1 The first spinor component of the \f$x\f$ vector. * @param xs2 The second spinor component of the \f$x\f$ vector. * @param xs3 The third spinor component of the \f$x\f$ vector. * @param xs4 The fourth spinor component of the \f$x\f$ vector. * @param ys1 The first spinor component of the \f$y\f$ vector. * @param ys2 The second spinor component of the \f$y\f$ vector. * @param ys3 The third spinor component of the \f$y\f$ vector. * @param ys4 The fourth spinor component of the \f$y\f$ vector. * @param zs1 The first spinor component of the \f$z\f$ vector. * @param zs2 The second spinor component of the \f$z\f$ vector. * @param zs3 The third spinor component of the \f$z\f$ vector. * @param zs4 The fourth spinor component of the \f$z\f$ vector. * @param ts1 The first spinor component of the \f$t\f$ vector. * @param ts2 The second spinor component of the \f$t\f$ vector. * @param ts3 The third spinor component of the \f$t\f$ vector. * @param ts4 The fourth spinor component of the \f$t\f$ vector. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex xs1, complex xs2, complex xs3, complex xs4, complex ys1, complex ys2, complex ys3, complex ys4, complex zs1, complex zs2, complex zs3, complex zs4, complex ts1, complex ts2, complex ts3, complex ts4) : WaveFunctionBase(p,part), _wf(xs1,xs2,xs3,xs4, ys1,ys2,ys3,ys4, zs1,zs2,zs3,zs4, ts1,ts2,ts3,ts4) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, const LorentzRSSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==4); } /** * Constructor, set the particle and the wavefunction. * @param p Particle * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorBarWaveFunction(const tPPtr & p, const LorentzRSSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==4); for (unsigned int i=0; i<4; ++i) for(unsigned int j=0; j<4; ++j) _wf(i,j)=wave(i,j)*UnitRemoval::InvSqrtE; } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1,2,3 as described above.) * @param dir The direction. */ RSSpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel,Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==4); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ RSSpinorBarWaveFunction(Lorentz5Momentum p,tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==4); } /** * Default constructor */ RSSpinorBarWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ RSSpinorBarWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * subscript operator for the wavefunction * Set components by index. */ complex operator ()(int i, int j) const { assert( i>=0 && i<=3 && j>=0 && j<=3 ); return _wf(i,j); } /** * return wavefunction as LorentzRSSpinorBar */ const LorentzRSSpinorBar & wave() const {return _wf;} + + /// Return wavefunction as LorentzRSSpinorBar + LorentzRSSpinorBar dimensionedWf() const { + LorentzRSSpinorBar temp(_wf.Type()); + for (unsigned int i=0; i<4; ++i) + for (unsigned int j=0; j<4; ++j) + temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; + return temp; + } /** * Get first spinor component for the x vector */ complex xs1() const {return _wf.xs1();} /** * Get second spinor component for the x vector */ complex xs2() const {return _wf.xs2();} /** * Get third spinor component for the x vector */ complex xs3() const {return _wf.xs3();} /** * Get fourth spinor component for the x vector */ complex xs4() const {return _wf.xs4();} /** * Get first spinor component for the y vector */ complex ys1() const {return _wf.ys1();} /** * Get second spinor component for the y vector */ complex ys2() const {return _wf.ys2();} /** * Get third spinor component for the y vector */ complex ys3() const {return _wf.ys3();} /** * Get fourth spinor component for the y vector */ complex ys4() const {return _wf.ys4();} /** * Get first spinor component for the z vector */ complex zs1() const {return _wf.zs1();} /** * Get second spinor component for the z vector */ complex zs2() const {return _wf.zs2();} /** * Get third spinor component for the z vector */ complex zs3() const {return _wf.zs3();} /** * Get fourth spinor component for the z vector */ complex zs4() const {return _wf.zs4();} /** * Get first spinor component for the t vector */ complex ts1() const {return _wf.ts1();} /** * Get second spinor component for the t vector */ complex ts2() const {return _wf.ts2();} /** * Get third spinor component for the t vector */ complex ts3() const {return _wf.ts3();} /** * Get fourth spinor component for the t vector */ complex ts4() const {return _wf.ts4();} //@} /** * reset functions */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1,2,3 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); + + /** + * Calculate the wavefunctions + */ + static void calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1,2,3 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * storage of the Lorentz RSSpinorBar */ LorentzRSSpinorBar _wf; - - /// Return wavefunction as LorentzRSSpinorBar - LorentzRSSpinorBar dimensionedWf() const { - LorentzRSSpinorBar temp(_wf.Type()); - for (unsigned int i=0; i<4; ++i) - for (unsigned int j=0; j<4; ++j) - temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; - return temp; - } }; } } #endif /* ThePEG_RSSpinorBarWaveFunction_H */ diff --git a/Helicity/WaveFunction/RSSpinorWaveFunction.cc b/Helicity/WaveFunction/RSSpinorWaveFunction.cc --- a/Helicity/WaveFunction/RSSpinorWaveFunction.cc +++ b/Helicity/WaveFunction/RSSpinorWaveFunction.cc @@ -1,269 +1,277 @@ // -*- 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" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace ThePEG; using namespace ThePEG::Helicity; // calculate the Wavefunction void RSSpinorWaveFunction::calculateWaveFunction(unsigned int ihel) { // if rest frame, make sure speical case is used Lorentz5Momentum ptemp=momentum(); double pr = ptemp.vect().mag2()/sqr(momentum().e()); if(pr<1e-20) { ptemp.setX(ZERO); ptemp.setY(ZERO); ptemp.setZ(ZERO); } assert(direction()!=intermediate); assert(ihel<=3); // only two valid helicities in massless case assert( mass()>ZERO || (ihel == 0 || ihel == 3 ) ); // new direct calculation _wf = LorentzRSSpinor(direction()==outgoing ? SpinorType::v : SpinorType::u); // compute the spinors std::array,2> spin; if(ihel!=0) spin[0] = HelicityFunctions::spinor(ptemp,1,direction()); if(ihel!=3) spin[1] = HelicityFunctions::spinor(ptemp,0,direction()); // compute the polarization vectors to construct the RS spinor std::array eps; if(ihel>=2) eps[0] = HelicityFunctions::polarizationVector(ptemp,2,direction(),vector_phase); else eps[2] = HelicityFunctions::polarizationVector(ptemp,0,direction(),vector_phase); if(mass()!=ZERO&&ihel!=0&&ihel!=3) eps[1] = HelicityFunctions::polarizationVector(ptemp,1,direction()); // 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(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = eps[0].x()*spin[0][iy]; _wf(1,iy) = eps[0].y()*spin[0][iy]; _wf(2,iy) = eps[0].z()*spin[0][iy]; _wf(3,iy) = eps[0].t()*spin[0][iy]; } } else if(ihel==2) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = or3*eps[0].x()*spin[1][iy]+tor3*eps[1].x()*spin[0][iy]; _wf(1,iy) = or3*eps[0].y()*spin[1][iy]+tor3*eps[1].y()*spin[0][iy]; _wf(2,iy) = or3*eps[0].z()*spin[1][iy]+tor3*eps[1].z()*spin[0][iy]; _wf(3,iy) = or3*eps[0].t()*spin[1][iy]+tor3*eps[1].t()*spin[0][iy]; } } else if(ihel==1) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = or3*eps[2].x()*spin[0][iy]+tor3*eps[1].x()*spin[1][iy]; _wf(1,iy) = or3*eps[2].y()*spin[0][iy]+tor3*eps[1].y()*spin[1][iy]; _wf(2,iy) = or3*eps[2].z()*spin[0][iy]+tor3*eps[1].z()*spin[1][iy]; _wf(3,iy) = or3*eps[2].t()*spin[0][iy]+tor3*eps[1].t()*spin[1][iy]; } } else if(ihel==0) { for(unsigned int iy=0;iy<4;++iy) { _wf(0,iy) = eps[2].x()*spin[1][iy]; _wf(1,iy) = eps[2].y()*spin[1][iy]; _wf(2,iy) = eps[2].z()*spin[1][iy]; _wf(3,iy) = eps[2].t()*spin[1][iy]; } } // this makes the phase choice the same as madgraph, useful for debugging only // Energy pt = ptemp.perp(); // double fact = direction()==incoming ? 1. : -1.; // Complex emphi(fact*ptemp.x()/pt,-fact*ptemp.y()/pt); // if(ihel==3) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) /= emphi; // } // else if(ihel==1) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) *= emphi; // } // else if(ihel==0) { // for(unsigned int ix=0;ix<4;++ix) // for(unsigned int iy=0;iy<4;++iy) // _wf(ix,iy) *= sqr(emphi); // } } 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; - } + calculateWaveFunctions(waves,particle->momentum(),particle->dataPtr(),dir); + } +} + +void RSSpinorWaveFunction:: +calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction dir) { + waves.resize(4); + RSSpinorWaveFunction wave(momentum,parton,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/RSSpinorWaveFunction.h b/Helicity/WaveFunction/RSSpinorWaveFunction.h --- a/Helicity/WaveFunction/RSSpinorWaveFunction.h +++ b/Helicity/WaveFunction/RSSpinorWaveFunction.h @@ -1,369 +1,376 @@ // -*- C++ -*- // // RSSpinorWaveFunction.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_RSSpinorWaveFunction_H #define ThePEG_RSSpinorWaveFunction_H // This is the declaration of the RSSpinorWaveFunction class. #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** \ingroup Helicity * * The RSSpinorWaveFunction class is designed to store the wavefunction * of a spin-3/2 particle in a form suitable for use in helicity amplitude * calculations of the matrix element using a similar philosophy to the * FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzRSSpinor class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and ParticleData for the fermion. * * This class also contains the code which does the actually calculation of the * spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$u\f$ spinor. * - outgoing calculates a \f$v\f$ spinor. * * The spinors are calculated using a Clebsch-Gordon decomposition in the rest-frame * for a massive particle and boosted to the lab-frame. For massless particles the * calculation is performed in the lab-frame (N.B. there are only two helicities * \f$\pm\frac32\f$ in this case.) * * N.B. In our convention 0 is the \f$-\frac32\f$ helicity state, * 1 is the \f$-\frac12\f$ helicity state, * 2 is the \f$+\frac12\f$ helicity state * 3 is the \f$+\frac32\f$ helicity state and * * @see WaveFunctionBase * @see LorentzRSSpinor * @see HelicityDefinitions * */ class RSSpinorWaveFunction: public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param xs1 The first spinor component of the \f$x\f$ vector. * @param xs2 The second spinor component of the \f$x\f$ vector. * @param xs3 The third spinor component of the \f$x\f$ vector. * @param xs4 The fourth spinor component of the \f$x\f$ vector. * @param ys1 The first spinor component of the \f$y\f$ vector. * @param ys2 The second spinor component of the \f$y\f$ vector. * @param ys3 The third spinor component of the \f$y\f$ vector. * @param ys4 The fourth spinor component of the \f$y\f$ vector. * @param zs1 The first spinor component of the \f$z\f$ vector. * @param zs2 The second spinor component of the \f$z\f$ vector. * @param zs3 The third spinor component of the \f$z\f$ vector. * @param zs4 The fourth spinor component of the \f$z\f$ vector. * @param ts1 The first spinor component of the \f$t\f$ vector. * @param ts2 The second spinor component of the \f$t\f$ vector. * @param ts3 The third spinor component of the \f$t\f$ vector. * @param ts4 The fourth spinor component of the \f$t\f$ vector. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex xs1, complex xs2, complex xs3, complex xs4, complex ys1, complex ys2, complex ys3, complex ys4, complex zs1, complex zs2, complex zs3, complex zs4, complex ts1, complex ts2, complex ts3, complex ts4) : WaveFunctionBase(p,part), _wf(xs1,xs2,xs3,xs4, ys1,ys2,ys3,ys4, zs1,zs2,zs3,zs4, ts1,ts2,ts3,ts4) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part, const LorentzRSSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==4); } /** * Constructor, set the momentum and the wavefunction. * @param p The Particle pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ RSSpinorWaveFunction(const tPPtr & p, const LorentzRSSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==4); for (unsigned int i=0; i<4; ++i) for(unsigned int j=0; j<4; ++j) _wf(i,j)=wave(i,j)*UnitRemoval::InvSqrtE; } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1,2,3 as described above.) * @param dir The direction. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel, Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==4); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ RSSpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==4); } /** * Default constructor */ RSSpinorWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ RSSpinorWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * subscript operator for the wavefunction * Set components by index. */ complex operator ()(int i, int j) const { assert( i>=0 && i<=3 && j>=0 && j<=3); return _wf(i,j); } /** * return wavefunction as LorentzRSSpinor */ const LorentzRSSpinor & wave() const {return _wf;} + /// Return wavefunction as LorentzRSSpinor + LorentzRSSpinor dimensionedWf() const { + LorentzRSSpinor temp(_wf.Type()); + for (unsigned int i=0; i<4; ++i) + for (unsigned int j=0; j<4; ++j) + temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; + return temp; + } + /** * Get first spinor component for the x vector */ complex xs1() const {return _wf.xs1();} /** * Get second spinor component for the x vector */ complex xs2() const {return _wf.xs2();} /** * Get third spinor component for the x vector */ complex xs3() const {return _wf.xs3();} /** * Get fourth spinor component for the x vector */ complex xs4() const {return _wf.xs4();} /** * Get first spinor component for the y vector */ complex ys1() const {return _wf.ys1();} /** * Get second spinor component for the y vector */ complex ys2() const {return _wf.ys2();} /** * Get third spinor component for the y vector */ complex ys3() const {return _wf.ys3();} /** * Get fourth spinor component for the y vector */ complex ys4() const {return _wf.ys4();} /** * Get first spinor component for the z vector */ complex zs1() const {return _wf.zs1();} /** * Get second spinor component for the z vector */ complex zs2() const {return _wf.zs2();} /** * Get third spinor component for the z vector */ complex zs3() const {return _wf.zs3();} /** * Get fourth spinor component for the z vector */ complex zs4() const {return _wf.zs4();} /** * Get first spinor component for the t vector */ complex ts1() const {return _wf.ts1();} /** * Get second spinor component for the t vector */ complex ts2() const {return _wf.ts2();} /** * Get third spinor component for the t vector */ complex ts3() const {return _wf.ts3();} /** * Get fourth spinor component for the t vector */ complex ts4() const {return _wf.ts4();} //@} /** * reset functions */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1,2,3 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ + static void calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction); + + /** + * Calculate the wavefunctions + */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1,2,3 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * storage of the Lorentz RSSpinor */ LorentzRSSpinor _wf; - - /// Return wavefunction as LorentzRSSpinor - LorentzRSSpinor dimensionedWf() const { - LorentzRSSpinor temp(_wf.Type()); - for (unsigned int i=0; i<4; ++i) - for (unsigned int j=0; j<4; ++j) - temp(i,j) = _wf(i,j)*UnitRemoval::SqrtE; - return temp; - } }; } } #endif /* ThePEG_RSSpinorWaveFunction_H */ diff --git a/Helicity/WaveFunction/SpinorBarWaveFunction.cc b/Helicity/WaveFunction/SpinorBarWaveFunction.cc --- a/Helicity/WaveFunction/SpinorBarWaveFunction.cc +++ b/Helicity/WaveFunction/SpinorBarWaveFunction.cc @@ -1,254 +1,260 @@ // -*- 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) { _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; - } + calculateWaveFunctions(waves,particle->momentum(),particle->dataPtr(),dir); } } +void SpinorBarWaveFunction::calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction dir) { + waves.resize(2); + SpinorBarWaveFunction wave(momentum,parton,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/SpinorBarWaveFunction.h b/Helicity/WaveFunction/SpinorBarWaveFunction.h --- a/Helicity/WaveFunction/SpinorBarWaveFunction.h +++ b/Helicity/WaveFunction/SpinorBarWaveFunction.h @@ -1,301 +1,308 @@ // -*- C++ -*- // // SpinorBarWaveFunction.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_SpinorBarWaveFunction_H #define ThePEG_SpinorBarWaveFunction_H // // This is the declaration of the SpinorBarWaveFunction class. #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** * Forward declaration of the SpinorWaveFunction class */ class SpinorWaveFunction; /** \ingroup Helicity * \author Peter Richardson * * The SpinorBarWaveFunction class is designed to store the wavefunction * of a barred spinor in a form suitable for use in helicity amplitude * calculations of the matrix element using a similar philosophy to the * FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzSpinorBar class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and ParticleData for the fermion. * * This class also contains the code which does the actually calculation * of the barred spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$\bar{v}\f$ spinor. * - outgoing calculates a \f$\bar{u}\f$ spinor. * * N.B. In our convention 0 is the \f$-\frac12\f$ helicity state and * 1 is the \f$+\frac12\f$ helicity state * * @see WaveFunctionBase * @see LorentzSpinorBar * @see HelicityDefinitions */ class SpinorBarWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param s1 The first component * @param s2 The second component * @param s3 The third component * @param s4 The fourth component */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, complex s1,complex s2, complex s3,complex s4) : WaveFunctionBase(p,part), _wf(s1,s2,s3,s4) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, const LorentzSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==2); } SpinorBarWaveFunction(const tPPtr & p, const LorentzSpinorBar & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(),p->dataPtr(),dir), _wf(wave.Type()) { assert(iSpin()==2); for (unsigned int i=0; i<4; ++i) _wf[i]=wave[i]*UnitRemoval::InvSqrtE; } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1 as described above.) * @param dir The direction. */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel,Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==2); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ SpinorBarWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==2); } /** * Default constructor. */ SpinorBarWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ SpinorBarWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Subscript operator for the wavefunction. */ complex operator ()(int i) const { assert(i>=0 &&i<=3); return _wf(i); } /** * Return wavefunction as LorentzSpinorBar. */ const LorentzSpinorBar & wave() const {return _wf;} /// Return wavefunction as LorentzSpinorBar LorentzSpinorBar dimensionedWave() const { return dimensionedWf(); } /** * Get the first spin component component. */ complex s1() const {return _wf.s1();} /** * Get the second spin component component. */ complex s2() const {return _wf.s2();} /** * Get the third spin component component. */ complex s3() const {return _wf.s3();} /** * Get the fourth spin component component. */ complex s4() const {return _wf.s4();} /** * Take the conjugate of the spinor \f$u_c=C\bar{u}^T\f$. This operation * transforms u-spinors to v-spinors and vice-versa and is required when * dealing with majorana particles. */ void conjugate(); /** * Return the barred spinor */ SpinorWaveFunction bar(); /** * Reset functions. */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1 as described above.) */ void calculateWaveFunction(unsigned int ihel); public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ + static void calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction); + + /** + * Calculate the wavefunctions + */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Storage of the Lorentz SpinorBar wavefunction. */ LorentzSpinorBar _wf; /// Return wavefunction as LorentzSpinorBar LorentzSpinorBar dimensionedWf() const { LorentzSpinorBar temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) temp(i) = _wf(i)*UnitRemoval::SqrtE; return temp; } }; } } #endif diff --git a/Helicity/WaveFunction/SpinorWaveFunction.cc b/Helicity/WaveFunction/SpinorWaveFunction.cc --- a/Helicity/WaveFunction/SpinorWaveFunction.cc +++ b/Helicity/WaveFunction/SpinorWaveFunction.cc @@ -1,250 +1,256 @@ // -*- 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" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace ThePEG; using namespace Helicity; // calculate the Wavefunction void SpinorWaveFunction::calculateWaveFunction(unsigned int ihel) { _wf = HelicityFunctions::spinor(momentum(),ihel,direction()); } 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; - } + calculateWaveFunctions(waves,particle->momentum(),particle->dataPtr(),dir); + } +} +void SpinorWaveFunction:: +calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction dir) { + SpinorWaveFunction wave(momentum,parton,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/Helicity/WaveFunction/SpinorWaveFunction.h b/Helicity/WaveFunction/SpinorWaveFunction.h --- a/Helicity/WaveFunction/SpinorWaveFunction.h +++ b/Helicity/WaveFunction/SpinorWaveFunction.h @@ -1,308 +1,315 @@ // -*- C++ -*- // // SpinorWaveFunction.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_SpinorWaveFunction_H #define ThePEG_SpinorWaveFunction_H // // This is the declaration of the SpinorWaveFunction class. // #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** * Forward declaration of the SpinorBarWaveFunction class */ class SpinorBarWaveFunction; /** \ingroup Helicity * \author Peter Richardson * * The SpinorWaveFunction class is designed to store the wavefunction * of a spinor in a form suitable for use in helicity amplitude calculations * of the matrix element using a similar philosophy to the FORTRAN HELAS code. * * In addition to storing the spinor using the LorentzSpinor class * it inherits from the WaveFunctionBase class to provide storage of * the momentum and getParticleData for the fermion. * * This class also contains the code which does the actually calculation * of the spinor for an external particle. * * When calculating the wavefunction the direction of the particle is used, * * \e i.e. * - incoming calculates a \f$u\f$ spinor. * - outgoing calculates a \f$v\f$ spinor. * * N.B. In our convention 0 is the \f$-\frac12\f$ helicity state and * 1 is the \f$+\frac12\f$ helicity state * * @see WaveFunctionBase * @see LorentzSpinor * @see HelicityDefinitions */ class SpinorWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and the components of the spinor. * @param p The momentum. * @param part The ParticleData pointer. * @param s1 The first component * @param s2 The second component * @param s3 The third component * @param s4 The fourth component */ SpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part,complex s1, complex s2,complex s3,complex s4) : WaveFunctionBase(p,part), _wf(s1,s2,s3,s4) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part, const LorentzSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==2); } /** * Constructor, set the momentum and the wavefunction. * @param p The particle * @param wave The wavefunction. * @param dir The direction of the particle */ SpinorWaveFunction(const tPPtr & p, const LorentzSpinor & wave, Direction dir=intermediate) : WaveFunctionBase(p->momentum(), p->dataPtr(), dir), _wf(wave.Type()) { assert(iSpin()==2); for (unsigned int i=0; i<4; ++i) _wf[i]=wave[i]*UnitRemoval::InvSqrtE; } /** * Constructor, set the momentum, helicity, direction. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1 as described above.) * @param dir The direction. */ SpinorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel, Direction dir) : WaveFunctionBase(p,part,dir) { assert(iSpin()==2); calculateWaveFunction(ihel); } /** * Constructor, set the momentum, direction, zero the * wavefunction. * @param p The momentum. * @param part The ParticleData pointer. * @param dir The direction. */ SpinorWaveFunction(const Lorentz5Momentum & p, tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==2); } /** * Default constructor. */ SpinorWaveFunction() : WaveFunctionBase(), _wf() {} /** * Special for spin correlations */ SpinorWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool=true) { calculateWaveFunctions(wave,part,dir); constructSpinInfo(wave,part,dir,time); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Subscript operator for the wavefunction. */ complex operator ()(int i) const { assert(i>=0 &&i<=3); return _wf(i); } /** * Return wavefunction as LorentzSpinor. */ const LorentzSpinor & wave() const {return _wf;} /** * Return wavefunction as LorentzSpinor. */ LorentzSpinor dimensionedWave() const {return dimensionedWf();} /** * Get the first spin component component. */ complex s1() const {return _wf.s1();} /** * Get the second spin component component. */ complex s2() const {return _wf.s2();} /** * Get the third spin component component. */ complex s3() const {return _wf.s3();} /** * Get the fourth spin component component. */ complex s4() const {return _wf.s4();} //@} /** * Take the conjugate of the spinor \f$u_c=C\bar{u}^T\f$. This operation * transforms u-spinors to v-spinors and vice-versa and is required when * dealing with majorana particles. */ void conjugate() { _wf=_wf.conjugate(); } /** * Return the barred spinor */ SpinorBarWaveFunction bar(); /** * Reset functions. */ //@{ /** * Reset the helicity (calculates the new spinor). * @param ihel The helicity (0,1 as described above.) */ void reset(unsigned int ihel) { calculateWaveFunction(ihel); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector > & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction); /** * Calculate the wavefunctions */ + static void calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton,Direction); + + /** + * Calculate the wavefunctions + */ static void calculateWaveFunctions(vector > & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector > & waves, tPPtr part,Direction dir, bool time); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time); private: /** * Calcuate the wavefunction. * @param ihel The helicity (0,1 as described above.) */ void calculateWaveFunction(unsigned int ihel); private: /** * Storage of the Lorentz Spinor. */ LorentzSpinor _wf; /// Return wavefunction as LorentzSpinor LorentzSpinor dimensionedWf() const { LorentzSpinor temp(_wf.Type()); for (unsigned int i=0; i<4; ++i) temp(i) = _wf(i)*UnitRemoval::SqrtE; return temp; } }; } } #endif diff --git a/Helicity/WaveFunction/VectorWaveFunction.cc b/Helicity/WaveFunction/VectorWaveFunction.cc --- a/Helicity/WaveFunction/VectorWaveFunction.cc +++ b/Helicity/WaveFunction/VectorWaveFunction.cc @@ -1,232 +1,240 @@ // -*- 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) { _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; - } + calculateWaveFunctions(waves,particle->momentum(),particle->dataPtr(), + dir,massless,phase); + } +} + +void VectorWaveFunction:: +calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, tcPDPtr parton, + Direction dir, bool massless, + VectorPhase phase) { + waves.resize(3); + VectorWaveFunction wave(momentum,parton,0,dir,phase); + for(unsigned int ix=0;ix<3;++ix) { + if(massless&&ix==1) { + waves[ix] = VectorWaveFunction(momentum,parton,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()); } } diff --git a/Helicity/WaveFunction/VectorWaveFunction.h b/Helicity/WaveFunction/VectorWaveFunction.h --- a/Helicity/WaveFunction/VectorWaveFunction.h +++ b/Helicity/WaveFunction/VectorWaveFunction.h @@ -1,255 +1,263 @@ // -*- C++ -*- // // VectorWaveFunction.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_VectorWaveFunction_H #define ThePEG_VectorWaveFunction_H // // This is the declaration of the VectorWaveFunction class. // #include "WaveFunctionBase.h" #include #include #include #include namespace ThePEG { namespace Helicity { /** \ingroup Helicity * * \author Peter Richardson * * The VectorWaveFunction class is designed to store the wavefunction * of a vector in a form suitable for use in helicity amplitude calculations * of the matrix element using a similar philosophy to the FORTRAN HELAS code. * * In addition to storing the vector using the LorentzPolarizationVector class * it inherits from the WaveFunctionBase class to provide storage of the * momentum and ParticleData for the vector boson. * * This class also contains the code which does the actually calculation of the * vector wavefunction. * * There are two choices available for the calculation of the wavefunction. * These are set using the VectorPhase enumeration which specifies a default choice. * The first choice, vector_phase, includes a phase factor \f$\exp(\pm i \phi)\f$ * for the \f$\pm\f$ helicity states while the second, vector_nophase, does not. * * N.B. In our convention 0 is the \f$-1\f$ helicity state and * 1 is the \f$0\f$ helicity state * 2 is the \f$+1\f$ helicity state * * @see WaveFunctionBase * @see LorentzPolarizationVector */ class VectorWaveFunction : public WaveFunctionBase { public: /** @name Standard constructors and destructors. */ //@{ /** * Constructor, set the momentum and Wavefunction, the direction can also * be specified. * @param p The momentum. * @param part The ParticleData pointer * @param wave The wavefunction, \e i.e. the polarization vector. * @param dir The direction of the particle. */ VectorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, const LorentzPolarizationVector & wave, Direction dir=intermediate) : WaveFunctionBase(p,part,dir), _wf(wave) { assert(iSpin()==3); } /** * Constructor, set the momentum and components of the wavefunction. * @param p The momentum. * @param part The ParticleData pointer * @param x The x component of the polarization vector * @param y The y component of the polarization vector * @param z The z component of the polarization vector * @param t The t component of the polarization vector */ VectorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part,const Complex & x, const Complex & y,const Complex & z, const Complex & t) : WaveFunctionBase(p,part), _wf(x,y,z,t) { assert(iSpin()==3); } /** * Constructor, set the momentum, helicity and direction, optionally the choice * of the phase. * @param p The momentum. * @param part The ParticleData pointer. * @param ihel The helicity (0,1,2 as described above.) * @param dir The direction. * @param phase The phase choice. */ VectorWaveFunction(const Lorentz5Momentum & p,tcPDPtr part, unsigned int ihel,Direction dir, VectorPhase phase=default_vector_phase) : WaveFunctionBase(p,part,dir) { assert(iSpin()==3); calculateWaveFunction(ihel,phase); } /** * Constructor, set the 5-momentum and direction, zero the wavefunction. * @param p The 5-momentum. * @param part The ParticleData pointer. * @param dir The direction. */ VectorWaveFunction(const Lorentz5Momentum &p, tcPDPtr part,Direction dir) : WaveFunctionBase(p,part,dir), _wf() { assert(iSpin()==3); } /** * Default constructor. */ VectorWaveFunction() {} /** * Special for spin correlations \todo make static? */ VectorWaveFunction(vector & wave, tPPtr part,Direction dir,bool time,bool massless, bool=true, VectorPhase phase=default_vector_phase) { calculateWaveFunctions(wave,part,dir,massless,phase); constructSpinInfo(wave,part,dir,time,massless); } //@} /** * Access to the wavefunction and its components. */ //@{ /** * Return wavefunction as polarization vector. */ const LorentzPolarizationVector & wave() const { return _wf;} /** * Get x component. */ Complex x() const {return _wf.x();} /** * Get y component. */ Complex y() const {return _wf.y();} /** * Get z component. */ Complex z() const {return _wf.z();} /** * Get t component. */ Complex t() const {return _wf.t();} /** * Reset functions. */ //@{ /** * Reset the helicity (recalculation the polarization vector). * @param ihel The new helicity (0,1,2 as described above.) * @param phase The phase choice. */ void reset(unsigned int ihel,VectorPhase phase=default_vector_phase) { calculateWaveFunction(ihel,phase); } //@} public: /** * Perform the Lorentz transformation of the wave function */ void transform(const LorentzRotation & r) { _wf.transform(r); transformMomentum(r); } public: /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction,bool massless, VectorPhase phase=default_vector_phase); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, tPPtr particle,Direction,bool massless, VectorPhase phase=default_vector_phase); /** * Calculate the wavefunctions */ + static void calculateWaveFunctions(vector & waves, + const Lorentz5Momentum & momentum, + tcPDPtr parton, Direction,bool massless, + VectorPhase phase=default_vector_phase); + + /** + * Calculate the wavefunctions + */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction,bool massless, VectorPhase phase=default_vector_phase); /** * Calculate the wavefunctions */ static void calculateWaveFunctions(vector & waves, RhoDMatrix & rho, tPPtr particle,Direction,bool massless, VectorPhase phase=default_vector_phase); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool massless); /** * Construct the SpinInfo object */ static void constructSpinInfo(const vector & waves, tPPtr part,Direction dir, bool time,bool massless); private: /** * Calculate the wavefunction * @param ihel The helicity (0,1,2 as described above.) * @param phase The phase choice. */ void calculateWaveFunction(unsigned int ihel, VectorPhase phase=default_vector_phase); private: /** * Storage of the wavefunction as a Lorentz Vector. */ LorentzPolarizationVector _wf; }; } } #endif