diff --git a/Decay/Baryon/Baryon1MesonDecayerBase.cc b/Decay/Baryon/Baryon1MesonDecayerBase.cc --- a/Decay/Baryon/Baryon1MesonDecayerBase.cc +++ b/Decay/Baryon/Baryon1MesonDecayerBase.cc @@ -1,961 +1,961 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the Baryon1MesonDecayerBase class. // #include "Baryon1MesonDecayerBase.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/PDT/DecayMode.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h" #include "ThePEG/Helicity/LorentzPolarizationVector.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" using namespace Herwig; using namespace ThePEG::Helicity; // The following static variable is needed for the type // description system in ThePEG. DescribeAbstractNoPIOClass describeHerwigBaryon1MesonDecayerBase("Herwig::Baryon1MesonDecayerBase", "HwBaryonDecay.so"); void Baryon1MesonDecayerBase::Init() { static ClassDocumentation documentation ("The Baryon1MesonDecayerBase class is the base class for" " the decays of the baryons to a baryon and a pseudoscalar or vector meson."); } void Baryon1MesonDecayerBase:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // for the decaying particle if(part.id()>0) { // incoming particle if(part.dataPtr()->iSpin()==PDT::Spin1Half) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&part),incoming,true); } else if(part.dataPtr()->iSpin()==PDT::Spin3Half) { RSSpinorWaveFunction:: constructSpinInfo(_inThreeHalf,const_ptr_cast(&part),incoming,true); } else assert(false); // outgoing fermion if(decay[0]->dataPtr()->iSpin()==PDT::Spin1Half) { SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else if(decay[0]->dataPtr()->iSpin()==PDT::Spin3Half) { RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else assert(false); } else { // incoming particle if(part.dataPtr()->iSpin()==PDT::Spin1Half) { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&part),incoming,true); } else if(part.dataPtr()->iSpin()==PDT::Spin3Half) { RSSpinorBarWaveFunction:: constructSpinInfo(_inThreeHalfBar,const_ptr_cast(&part),incoming,true); } else assert(false); // outgoing fermion if(decay[0]->dataPtr()->iSpin()==PDT::Spin1Half) { SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } else if(decay[0]->dataPtr()->iSpin()==PDT::Spin3Half) { RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } else assert(false); } // outgoing meson if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) { ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); } else if(decay[1]->dataPtr()->iSpin()==PDT::Spin1) { VectorWaveFunction::constructSpinInfo(_inVec,decay[1],outgoing,true, decay[1]->id()==ParticleID::gamma); } else assert(false); } // return the matrix element squared for a given mode and phase-space channel // (inherited from DecayIntegrator and implemented here) double Baryon1MesonDecayerBase::me2(const int ichan,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { double me(0.); // decide which matrix element we are doing // incoming spin-1/2 particle if(part.dataPtr()->iSpin()==2) { // decay to spin-1/2 particle if(outgoing[0]->iSpin()==2) { // scalar meson if(outgoing[1]->iSpin()==1) me=halfHalfScalar(ichan,part,outgoing,momenta,meopt); // vector meson else if(outgoing[1]->iSpin()==3) me=halfHalfVector(ichan,part,outgoing,momenta,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // decay to spin-3/2 particle else if(outgoing[0]->iSpin()==4) { // scalar meson if(outgoing[1]->iSpin()==1) me=halfThreeHalfScalar(ichan,part,outgoing,momenta,meopt); // vector meson else if(outgoing[1]->iSpin()==3) me=halfThreeHalfVector(ichan,part,outgoing,momenta,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // unknown else throw DecayIntegratorError() << "Unknown outgoing baryon spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // incoming spin-3/2 particle else if(part.dataPtr()->iSpin()==4) { // decay to spin-1/2 particle if(outgoing[0]->iSpin()==2) { // scalar meson if(outgoing[1]->iSpin()==1) me=threeHalfHalfScalar(ichan,part,outgoing,momenta,meopt); // vector meson else if(outgoing[1]->iSpin()==3) me=threeHalfHalfVector(ichan,part,outgoing,momenta,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // decay to spin-3/2 particle else if(outgoing[0]->iSpin()==4) { // scalar meson if(outgoing[1]->iSpin()==1) me=threeHalfThreeHalfScalar(ichan,part,outgoing,momenta,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // unknown else throw DecayIntegratorError() << "Unknown outgoing baryon spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // unknown else throw DecayIntegratorError() << "Unknown incoming spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; return me; } // matrix element for the decay of a spin-1/2 fermion to a spin-1/2 fermion and // a pseudoscalar meson double Baryon1MesonDecayerBase:: halfHalfScalar(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // matrix element } // spinors for the decay product if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalfBar[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalf[ix] = HelicityFunctions::dimensionedSpinor (-momenta[0],ix,Helicity::outgoing); } // get the couplings Complex A,B; halfHalfScalarCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(),A,B); Complex left,right,meout; // coupling for an incoming particle if(part.id()>0) { left = (A-B); right = (A+B); } // coupling for an incoming antiparticle else { left = conj(A+B); right = conj(A-B); } // calculate the matrix element vector ispin(3,0); unsigned int ix,iy; // Complex output(0.); for(ix=0;ix<2;++ix) { for(iy=0;iy<2;++iy) { if(outgoing[0]->id()>0){ispin[0]=iy;ispin[1]=ix;} else{ispin[0]=ix;ispin[1]=iy;} (*ME())(ispin)=Complex(_inHalf[iy].generalScalar(_inHalfBar[ix],left,right)/part.mass()); } } // test of the matrix element // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // Complex h1(2.*Qp*A/part.mass()),h2(-2.*Qm*B/part.mass()); // cout << "testing 1/2->1/2 0 " // << 0.5*output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)) << " " // << 0.5*(h1*conj(h1)+h2*conj(h2))/output << endl; // cout << "testing alpha " << // (norm(0.5*(h1+h2))-norm(0.5*(h1-h2)))/ // (norm(0.5*(h1+h2))+norm(0.5*(h1-h2))) << "\n"; // store the matrix element return (ME()->contract(_rho)).real(); } // matrix element for the decay of a spin-1/2 fermion to a spin-1/2 fermion and // a vector meson double Baryon1MesonDecayerBase:: halfHalfVector(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1))); // check if the outgoing meson is really a photon bool photon=outgoing[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // matrix element } // spinors for the decay product if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalfBar[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalf[ix] = HelicityFunctions::dimensionedSpinor (-momenta[0],ix,Helicity::outgoing); } _inVec.resize(3); for(unsigned int ix=0;ix<3;++ix) { if(photon && ix==1) continue; _inVec[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing); } // get the couplings Complex A1,A2,B1,B2; halfHalfVectorCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(), A1,A2,B1,B2); Complex lS,rS,lV,rV; complex scalar; // couplings for an incoming particle if(part.id()>0) { lS = (A2-B2); rS = (A2+B2); lV = (A1-B1); rV = (A1+B1); } else { lS = -conj(A2+B2); rS = -conj(A2-B2); lV = conj(A1-B1); rV = conj(A1+B1); } // calculate the matrix element // decide which type of mode to do Energy msum(part.mass()+momenta[0].mass()); vector ispin(3); LorentzVector > svec; Complex prod; // Complex output(0.); unsigned int ix,iy; for(ix=0;ix<2;++ix) { for(iy=0;iy<2;++iy) { // scalar like piece scalar = _inHalf[iy].generalScalar(_inHalfBar[ix],lS,rS); // vector like piece svec = _inHalf[iy].generalCurrent(_inHalfBar[ix],lV,rV); if(outgoing[0]->id()>0) { ispin[0] = iy; ispin[1] = ix; } else { ispin[0] = ix; ispin[1] = iy; } for(ispin[2]=0;ispin[2]<3;++ispin[2]) { ispin[2]=ispin[2]; prod=_inVec[ispin[2]].dot(part.momentum())/msum; (*ME())(ispin)=(svec.dot(_inVec[ispin[2]])+prod*scalar)/part.mass(); // output += norm((*ME())(ispin)); } } } // test of the matrix element // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // double r2(sqrt(2.)); // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // if(m3!=ZERO) { // Complex h1(2.*r2*Qp*B1/part.mass()),h2(-2.*r2*Qm*A1/part.mass()), // h3(2./m3*(Qp*(m1-m2)*B1-Qm*m1*B2*pcm/(m1+m2))/part.mass()), // h4(2./m3*(Qm*(m1+m2)*A1+Qp*m1*A2*pcm/(m1+m2))/part.mass()); // generator()->log() << "testing 1/2->1/2 1 " // << 0.5*output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+h4*conj(h4)) << " " // << 0.50*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+h4*conj(h4))/output // << "\n"; // generator()->log() << "alpha = " << 2.*(norm(h3)+norm(h4))/(norm(h1)+norm(h2))-1. // << "\n"; // generator()->log() << "masses " << m1/GeV << " " << m2/GeV << " " << m3/GeV << "\n"; // } // return the answer return (ME()->contract(_rho)).real(); } // matrix element for the decay of a spin-1/2 fermion to a spin-3/2 fermion and // a scalar meson double Baryon1MesonDecayerBase:: halfThreeHalfScalar(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin3Half,PDT::Spin0))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // matrix element } // spinors for the decay product LorentzPolarizationVector in=UnitRemoval::InvE*part.momentum(); if(part.id()>0) { RSSpinorBarWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalfBar.resize(4); _inHalfBar.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalfBar[ihel]=swave.dimensionedWf(); _inHalfBar[ihel] = _inThreeHalfBar[ihel].dot(in); } } else { RSSpinorWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalf.resize(4); _inHalf.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalf[ihel]=swave.dimensionedWf(); _inHalf[ihel] = _inThreeHalf[ihel].dot(in); } } // get the couplings Complex A,B,left,right; Energy msum(part.mass()+momenta[0].mass()); halfThreeHalfScalarCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(), A,B); // incoming particle if(part.id()>0) { left=(A-B); right=(A+B); } // incoming anti-particle else { left=conj(A+B); right=conj(A-B); } vector ispin(3,0); for(unsigned ixa=0;ixa<2;++ixa) { for(unsigned int iya=0;iya<4;++iya) { unsigned int ix(iya),iy(ixa); if(outgoing[0]->id()<0) swap(ix,iy); ispin[0]=ixa; ispin[1]=iya; complex value = _inHalf[iy].generalScalar(_inHalfBar[ix],left,right) *UnitRemoval::E/part.mass()/msum; (*ME())(ispin) = value; } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // double r23(sqrt(2./3.)); // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // complex h1(-2.*r23*pcm*m1/m2*Qm*B/(m1+m2)),h2( 2.*r23*pcm*m1/m2*Qp*A/(m1+m2)); // cout << "testing 1/2->3/2 0 " << part.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass())/output << endl; // return the answer return output; } // matrix element for the decay of a spin-1/2 fermion to a spin-3/2 fermion and // a vector meson double Baryon1MesonDecayerBase:: halfThreeHalfVector(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1Half,PDT::Spin3Half,PDT::Spin1))); // check if the outgoing meson is really a photon bool photon=outgoing[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&part), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&part), incoming); // matrix element } LorentzPolarizationVector in=UnitRemoval::InvE*part.momentum(); // wavefunctions for outgoing fermion if(part.id()>0) { RSSpinorBarWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalfBar.resize(4); _inHalfBar.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalfBar[ihel]=swave.dimensionedWf(); _inHalfBar[ihel] = _inThreeHalfBar[ihel].dot(in); } } else { RSSpinorWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalf.resize(4); _inHalf.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalf[ihel]=swave.dimensionedWf(); _inHalf[ihel] = _inThreeHalf[ihel].dot(in); } } ME()->zero(); _inVec.resize(3); for(unsigned int ix=0;ix<3;++ix) { if(photon && ix==1) continue; _inVec[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing); } // get the couplings Complex A1,A2,A3,B1,B2,B3; halfThreeHalfVectorCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(), A1,A2,A3,B1,B2,B3); Energy msum(part.mass()+momenta[0].mass()); Complex lS,rS,lV,rV,left,right; // incoming particle if(part.id()>0) { lS=(A3-B3);rS=(A3+B3); lV=(A2-B2);rV=(A2+B2); left=(A1-B1);right=(A1+B1); } // incoming anti-particle else { lS=conj(A3+B3);rS=conj(A3-B3); lV=-conj(A2-B2);rV=-conj(A2+B2); left=conj(A1+B1);right=conj(A1-B1); } // compute the matrix element vector ispin(3); LorentzVector > svec; Complex prod; complex scalar; LorentzSpinor stemp; LorentzSpinorBar sbtemp; for(unsigned iya=0;iya<4;++iya) { ispin[1]=iya; // piece where the vector-spinor is dotted with the momentum of the // incoming fermion for(unsigned ixa=0;ixa<2;++ixa) { unsigned int ix(iya),iy(ixa); if(outgoing[0]->id()<0) swap(ix,iy); scalar = _inHalf[iy].generalScalar (_inHalfBar[ix],lS,rS); svec = _inHalf[iy].generalCurrent(_inHalfBar[ix],lV,rV); ispin[0]=ixa; for(unsigned int iz=0;iz<3;++iz) { ispin[2]=iz; prod=_inVec[iz].dot(part.momentum())/msum; (*ME())(ispin) += (svec.dot(_inVec[iz])+prod*scalar)* UnitRemoval::E/msum/part.mass(); } } // the piece where the vector spinor is dotted with the polarization vector for(unsigned int iz=0;iz<3;++iz) { ispin[2]=iz; if(outgoing[0]->id()>0) sbtemp = _inThreeHalfBar[iya].dot(_inVec[iz]); else stemp = _inThreeHalf[iya].dot(_inVec[iz]); for(unsigned int ixa=0;ixa<2;++ixa) { ispin[0]=ixa; if(outgoing[0]->id()>0) stemp = _inHalf[ixa]; else (*ME())(ispin) += Complex(stemp.generalScalar(sbtemp,left,right)/part.mass()); } } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy2 m12(m1*m1),m22(m2*m2),m32(m3*m3); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // double r2(sqrt(2.)),r3(sqrt(3.)); // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // complex h1(-2.*Qp*A1),h2(2.*Qm*B1); // complex h3(-2./r3*Qp*(A1-Qm*Qm/m2*A2/msum)); // complex h4( 2./r3*Qm*(B1-Qp*Qp/m2*B2/msum)); // complex h5(-2.*r2/r3/m2/m3*Qp*(0.5*(m12-m22-m32)*A1+0.5*Qm*Qm*(m1+m2)*A2/msum // +m12*pcm*pcm*A3/msum/msum)); // complex h6( 2.*r2/r3/m2/m3*Qm*(0.5*(m12-m22-m32)*B1-0.5*Qp*Qp*(m1-m2)*B2/msum // +m12*pcm*pcm*B3/msum/msum)); // cout << "testing 1/2->3/2 1 " << part.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(part.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(part.mass())/output << endl; // return the answer return output; } // matrix element for the decay of a spin-3/2 fermion to a spin-1/2 fermion and // a scalar meson double Baryon1MesonDecayerBase:: threeHalfHalfScalar(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin3Half,PDT::Spin1Half,PDT::Spin0))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&part), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&part), incoming); } } // spinors for the decay product if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalfBar[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalf[ix] = HelicityFunctions::dimensionedSpinor (-momenta[0],ix,Helicity::outgoing); } LorentzPolarizationVector out=UnitRemoval::InvE*momenta[0]; if(part.id()>0) { _inHalf.resize(_inThreeHalf.size()); for(unsigned int ix=0;ix<_inThreeHalf.size();++ix) _inHalf[ix] = _inThreeHalf[ix].dot(out); } else { _inHalfBar.resize(_inThreeHalfBar.size()); for(unsigned int ix=0;ix<_inThreeHalfBar.size();++ix) _inHalfBar[ix] = _inThreeHalfBar[ix].dot(out); } // get the couplings Complex A,B; Energy msum=part.mass()+momenta[0].mass(); threeHalfHalfScalarCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(), A,B); Complex left,right; // incoming particle if(part.id()>0) { left=(A-B); right=(A+B); } // incoming anti-particle else { left=conj(A+B); right=conj(A-B); } // compute the matrix element vector ispin(3,0); for(unsigned ixa=0;ixa<2;++ixa) { for(unsigned iya=0;iya<4;++iya) { unsigned int iy=iya,ix=ixa; if(outgoing[0]->id()<0) swap(ix,iy); ispin[0]=iya; ispin[1]=ixa; (*ME())(ispin) = Complex(_inHalf[iy].generalScalar(_inHalfBar[ix],left,right)* UnitRemoval::E/msum/part.mass()); } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // double r23(sqrt(2./3.)); // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // complex h1(-2.*r23*pcm*Qm*B/(m1+m2)),h2( 2.*r23*pcm*Qp*A/(m1+m2)); - // cout << "testing 3/2->1/2 0 " << part.id() << " " - // << output << " " - // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass()) << " " - // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass())/output << endl; + // generator()->log() << "testing 3/2->1/2 0 " << part.id() << " " + // << output << " " + // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass()) << " " + // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(part.mass())/output << endl; // return the answer return output; } // matrix element for the decay of a spin-3/2 fermion to a spin-3/2 fermion and // a scalar meson double Baryon1MesonDecayerBase::threeHalfThreeHalfScalar(const int, const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin3Half,PDT::Spin3Half,PDT::Spin0))); // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&part), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&part), incoming); } // matrix element } // spinors for the decay product // wavefunctions for outgoing fermion if(part.id()>0) { RSSpinorBarWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalfBar.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalfBar[ihel]=swave.dimensionedWf(); } } else { RSSpinorWaveFunction swave(momenta[0],outgoing[0],Helicity::outgoing); _inThreeHalf.resize(4); for(unsigned int ihel=0;ihel<4;++ihel) { swave.reset(ihel); _inThreeHalf[ihel]=swave.dimensionedWf(); } } LorentzPolarizationVector in = UnitRemoval::InvE*part.momentum(); _inHalf.resize(_inThreeHalf.size()); _inHalfBar.resize(_inThreeHalfBar.size()); for(unsigned int ix=0;ix<_inThreeHalf.size();++ix) { _inHalf[ix] = _inThreeHalf[ix].dot(in); _inHalfBar[ix] = _inThreeHalfBar[ix].dot(in); } // get the couplings Complex A1,B1,A2,B2; Energy msum(part.mass()+momenta[0].mass()); threeHalfThreeHalfScalarCoupling(imode(),part.mass(),momenta[0].mass(), momenta[1].mass(),A1,A2,B1,B2); Complex left1,right1,left2,right2; // incoming particle if(part.id()>0) { left1=(A1-B1); right1=(A1+B1); left2=(A2-B2); right2=(A2+B2); } // incoming anti-particle else { left1=(A1+B1); right1=(A1-B1); left2=(A2+B2); right2=(A2-B2); } // compute the matrix element vector ispin(3,0); for(unsigned ixa=0;ixa<4;++ixa) { for(unsigned iya=0;iya<4;++iya) { unsigned int iy=iya,ix=ixa; if(outgoing[0]->id()<0) swap(ix,iy); ispin[0]=iya; ispin[1]=ixa; (*ME())(ispin)=Complex((_inThreeHalf[iy].generalScalar(_inThreeHalfBar[ix],left1,right1) +_inHalf[iy].generalScalar( _inHalfBar[ix],left2,right2) *UnitRemoval::E2/sqr(msum))/part.mass()); } } // return the answer return (ME()->contract(_rho)).real(); } // matrix element for the decay of a spin-3/2 fermion to a spin-1/2 fermion and // a vector meson double Baryon1MesonDecayerBase:: threeHalfHalfVector(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin3Half,PDT::Spin1Half,PDT::Spin1))); // check if the outgoing meson is really a photon bool photon=outgoing[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(part.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&part), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&part), incoming); } } if(part.id()>0) { _inHalfBar.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalfBar[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ix,Helicity::outgoing); } else { _inHalf.resize(2); for(unsigned int ix=0;ix<2;++ix) _inHalf[ix] = HelicityFunctions::dimensionedSpinor (-momenta[0],ix,Helicity::outgoing); } LorentzPolarizationVector out=UnitRemoval::InvE*momenta[0]; if(part.id()>0) { _inHalf.resize(_inThreeHalf.size()); for(unsigned int ix=0;ix<_inThreeHalf.size();++ix) _inHalf[ix] = _inThreeHalf[ix].dot(out); } else { _inHalfBar.resize(_inThreeHalfBar.size()); for(unsigned int ix=0;ix<_inThreeHalfBar.size();++ix) _inHalfBar[ix] = _inThreeHalfBar[ix].dot(out); } // wavefunctions for outgoing fermion ME()->zero(); _inVec.resize(3); for(unsigned int ix=0;ix<3;++ix) { if(photon && ix==1) continue; _inVec[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing); } // get the couplings Complex A1,A2,A3,B1,B2,B3,prod,meout; threeHalfHalfVectorCoupling(imode(),part.mass(),momenta[0].mass(),momenta[1].mass(), A1,A2,A3,B1,B2,B3); Energy msum(part.mass()+momenta[0].mass()); Complex lS,rS,lV,rV,left,right; // incoming particle if(part.id()>0) { lS=(A3-B3);rS=(A3+B3); lV=(A2-B2);rV=(A2+B2); left=(A1-B1);right=(A1+B1); } // incoming anti-particle else { lS=conj(A3+B3);rS=conj(A3-B3); lV=-conj(A2-B2);rV=-conj(A2+B2); left=conj(A1+B1);right=conj(A1-B1); } // compute the matrix element vector ispin(3); LorentzVector > svec; LorentzSpinor stemp; LorentzSpinorBar sbtemp; complex scalar; for(unsigned iya=0;iya<4;++iya) { ispin[0]=iya; for(unsigned ixa=0;ixa<2;++ixa) { unsigned int iy=iya,ix=ixa; if(outgoing[0]->id()<0) swap(ix,iy); scalar = _inHalf[iy].generalScalar( _inHalfBar[ix],lS,rS); svec = _inHalf[iy].generalCurrent(_inHalfBar[ix],lV,rV); ispin[1]=ixa; for(unsigned int iz=0;iz<3;++iz) { ispin[2]=iz; prod=_inVec[iz].dot(momenta[0])/msum; (*ME())(ispin) += (svec.dot(_inVec[iz])+prod*scalar)* UnitRemoval::E/msum/part.mass(); } } // the piece where the vector spinor is dotted with the polarization vector for(unsigned iz=0;iz<3;++iz) { ispin[2]=iz; if(outgoing[0]->id()>0) stemp = _inThreeHalf[iya].dot(_inVec[iz]); else sbtemp = _inThreeHalfBar[iya].dot(_inVec[iz]); for(unsigned int ixa=0;ixa<2;++ixa) { ispin[1]=ixa; if(outgoing[0]->id()>0) sbtemp = _inHalfBar[ixa]; else stemp = _inHalf[ixa]; (*ME())(ispin) += Complex(stemp.generalScalar(sbtemp,left,right)/part.mass()); } } } double output = (ME()->contract(_rho)).real(); // testing code // Energy m1(part.mass()),m2(momenta[0].mass()),m3(momenta[1].mass()); // Energy2 m12(m1*m1),m22(m2*m2),m32(m3*m3); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // double r2(sqrt(2.)),r3(sqrt(3.)); // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // complex h1(-2.*Qp*A1),h2(2.*Qm*B1); // complex h3(-2./r3*Qp*(A1-Qm*Qm/m1*A2/msum)); // complex h4( 2./r3*Qm*(B1-Qp*Qp/m1*B2/msum)); // complex h5(ZERO),h6(ZERO); // if(m3!=ZERO) { // h5 = (-2.*r2/r3/m1/m3*Qp*(0.5*(m22-m12-m32)*A1+0.5*Qm*Qm*(m1+m2)*A2/msum // +m12*pcm*pcm*A3/msum/msum)); // h6 = ( 2.*r2/r3/m1/m3*Qm*(0.5*(m22-m12-m32)*B1-0.5*Qp*Qp*(m2-m1)*B2/msum // +m22*pcm*pcm*B3/msum/msum)); // } // cout << "testing 3/2->1/2 1 " << part.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(part.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(part.mass())/output << endl; // return the answer return output; } void Baryon1MesonDecayerBase::halfHalfScalarCoupling(int,Energy,Energy,Energy, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::halfHalfScalarCoupling()" << " called from base class this must be implemented" << " in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::halfHalfVectorCoupling(int,Energy,Energy,Energy, Complex&,Complex&, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::halfHalfVectorCoupling()" << " called from base class this must be implemented " << "in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::halfThreeHalfScalarCoupling(int,Energy,Energy,Energy, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::halfThreeHalfScalarCoupling" << "() called from base class this must be implemented" << " in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::halfThreeHalfVectorCoupling(int,Energy,Energy,Energy, Complex&,Complex&, Complex&,Complex&, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::halfThreeHalfVectorCoupling" << "() called from base class this must be implemented " << "in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::threeHalfHalfScalarCoupling(int,Energy,Energy,Energy, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::threeHalfHalfScalarCoupling" << "() called from base class this must be implemented" << " in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::threeHalfHalfVectorCoupling(int,Energy,Energy,Energy, Complex&,Complex&, Complex&,Complex&, Complex&,Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::threeHalfHalfVectorCoupling" << "() called from base class this must be implemented " << "in the inheriting class" << Exception::abortnow; } void Baryon1MesonDecayerBase::threeHalfThreeHalfScalarCoupling(int,Energy,Energy, Energy,Complex&, Complex&,Complex&, Complex&) const { throw DecayIntegratorError() << "Baryon1MesonDecayerBase::threeHalfThreeHalfScalar" << "Coupling() called from base class this must be " << "implemented in the inheriting class" << Exception::abortnow; } bool Baryon1MesonDecayerBase::twoBodyMEcode(const DecayMode & dm,int & mecode, double & coupling) const { coupling=1.; unsigned int inspin(dm.parent()->iSpin()),outspin,outmes; ParticleMSet::const_iterator pit(dm.products().begin()); bool order; if((**pit).iSpin()%2==0) { order=true; outspin=(**pit).iSpin(); ++pit;outmes=(**pit).iSpin(); } else { order=false; outmes=(**pit).iSpin();++pit; outspin=(**pit).iSpin(); } mecode=-1; if(inspin==2) { if(outspin==2){if(outmes==1){mecode=101;}else{mecode=102;}} else if(outspin==4){if(outmes==1){mecode=103;}else{mecode=104;}} } else if(inspin==4) { if(outspin==2){if(outmes==1){mecode=105;}else{mecode=106;}} else if(outspin==4){if(outmes==1){mecode=107;}else{mecode=108;}} } return order; } void Baryon1MesonDecayerBase::dataBaseOutput(ofstream & os,bool header) const { DecayIntegrator::dataBaseOutput(os,header); } diff --git a/Decay/Partonic/WeakPartonicDecayer.cc b/Decay/Partonic/WeakPartonicDecayer.cc --- a/Decay/Partonic/WeakPartonicDecayer.cc +++ b/Decay/Partonic/WeakPartonicDecayer.cc @@ -1,611 +1,637 @@ // -*- C++ -*- // // WeakPartonicDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2002-2017 The Herwig Collaboration // // Herwig 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 WeakPartonicDecayer class. // #include "WeakPartonicDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "Herwig/Utilities/Kinematics.h" #include "ThePEG/PDT/EnumParticles.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Interface/Reference.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/PDT/ConstituentParticleData.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" using namespace Herwig; using namespace ThePEG::Helicity; WeakPartonicDecayer::WeakPartonicDecayer() : MECode(0), _radprob(0.0), _maxtry(300), _threemax(3.), _fourmax(3.) {} IBPtr WeakPartonicDecayer::clone() const { return new_ptr(*this); } IBPtr WeakPartonicDecayer::fullclone() const { return new_ptr(*this); } bool WeakPartonicDecayer::accept(tcPDPtr parent, const tPDVector & prod) const { // check we can find the flavours of the quarks in the decaying meson long id = parent->id(); int flav1, flav2; if((id / 1000)%10) { flav1 = (id/1000)%10; flav2 = (id/10)%100; } else { flav1 = id/100; flav2 = (id/10)%10; } if(!flav1 || !flav2) return false; // if two decay products one must be in triplet and one antitriplet if(prod.size()==2) { if((prod[0]->iColour()==PDT::Colour3&&prod[1]->iColour()==PDT::Colour3bar)|| (prod[0]->iColour()==PDT::Colour3bar&&prod[1]->iColour()==PDT::Colour3)) return true; } else if(prod.size()==3) { if(((prod[0]->iColour()==PDT::Colour3 &&prod[2]->iColour()==PDT::Colour3bar)|| (prod[0]->iColour()==PDT::Colour3bar&&prod[2]->iColour()==PDT::Colour3)) &&prod[1]->iColour()==PDT::Colour8) return true; } else if(prod.size()==4) { // first two particles should be leptons or q qbar if((prod[0]->id()>=11&&prod[0]->id()<=16&&prod[1]->id()<=-11&&prod[1]->id()>=-16)|| (prod[1]->id()>=11&&prod[1]->id()<=16&&prod[0]->id()<=-11&&prod[0]->id()>=-16)|| (prod[0]->iColour()==PDT::Colour3 &&prod[1]->iColour()==PDT::Colour3bar )|| (prod[1]->iColour()==PDT::Colour3 &&prod[0]->iColour()==PDT::Colour3bar )) { // third particle quark and fourth colour anti-triplet or // thrid particle antiquark and fourth colour triplet if((prod[2]->iColour()==PDT::Colour3bar&&prod[3]->iColour()==PDT::Colour3 )|| (prod[2]->iColour()==PDT::Colour3 &&prod[3]->iColour()==PDT::Colour3bar)) return true; } else return false; } return false; } ParticleVector WeakPartonicDecayer::decay(const Particle & parent, const tPDVector & children) const { static tcPDPtr gluon=getParticleData(ParticleID::g); // make the particles ParticleVector partons; for(unsigned int ix=0;ixproduceParticle()); // these products have the mass but should have constituent mass partons[ix]->set5Momentum(Lorentz5Momentum(children[ix]->constituentMass())); } // 2-body decays if(partons.size()==2) { // no gluon if not select based on probability or if three body not allowed if(UseRandom::rnd()>_radprob|| parent.mass()constituentMass()+partons[0]->mass()+partons[1]->mass()) { double ctheta,phi; Lorentz5Momentum pout[2]; for(unsigned int ix=0;ix<2;++ix) pout[ix].setMass(partons[ix]->mass()); Kinematics::generateAngles(ctheta,phi); Kinematics::twoBodyDecay(parent.momentum(),pout[0].mass(),pout[1].mass(), ctheta,phi,pout[0],pout[1]); for(unsigned int ix=0; ix<2;++ix) partons[ix]->setMomentum(pout[ix]); if(partons[0]->dataPtr()->iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[1]); } else { partons[0]-> colourNeighbour(partons[1]); } } else { Lorentz5Momentum pout[3]; for(unsigned int ix=0;ix<2;++ix) pout[ix].setMass(partons[ix]->mass()); // add gluon partons.push_back(gluon->produceParticle()); partons.back()->set5Momentum(gluon->constituentMass()); // momentum of gluon pout[2] = Lorentz5Momentum(gluon->constituentMass()); Kinematics::threeBodyDecay(parent.momentum(),pout[1],pout[0],pout[2]); for(unsigned int ix=0; ix<3;++ix) partons[ix]->setMomentum(pout[ix]); if(partons[0]->dataPtr()->iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[2]); partons[1]-> colourNeighbour(partons[2]); } else { partons[0]-> colourNeighbour(partons[2]); partons[1]->antiColourNeighbour(partons[2]); } } } // 3-body decays else if(partons.size()==3) { // set masses of products Lorentz5Momentum pout[3],pin(parent.momentum()); for(unsigned int ix=0;ix<3;++ix) pout[ix].setMass(partons[ix]->mass()); double xs(partons[2]->mass()/pin.e()),xb(1.-xs); pout[2]=xs*pin; // Get the particle quark that is decaying long idQ, idSpec; idSpec = partons[2]->id(); idQ = (parent.id()/1000)%10; if(!idQ) idQ = (parent.id()/100)%10; // Now the odd case of a B_c where the c decays, not the b if(idSpec == idQ) idQ = (parent.id()/10)%10; // momentum of the decaying quark PPtr inter = getParticleData(idQ)->produceParticle(parent.momentum()*xb); // two body decay of heavy quark double ctheta,phi; Kinematics::generateAngles(ctheta,phi); Kinematics::twoBodyDecay(inter->momentum(),pout[0].mass(),pout[1].mass(), ctheta,phi,pout[0],pout[1]); // set the momenta of the decay products for(unsigned int ix=0; ix<3;++ix) partons[ix]->setMomentum(pout[ix]); // make the colour connections // quark first if(partons[0]->data().iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[1]); partons[1]->colourNeighbour(partons[0]); partons[1]->antiColourNeighbour(partons[2]); partons[2]->colourNeighbour(partons[1]); } // antiquark first else { partons[0]->colourNeighbour(partons[1]); partons[1]->antiColourNeighbour(partons[0]); partons[1]->colourNeighbour(partons[2]); partons[2]->antiColourNeighbour(partons[1]); } } // 4-body decays else if(partons.size()==4) { // swap 0 and 1 if needed if(partons[1]->dataPtr()->iColour()!=PDT::Colour0&& partons[1]->dataPtr()->iColour()!=partons[2]->dataPtr()->iColour()) swap(partons[0],partons[1]); // get the momenta of the decaying quark and the spectator Lorentz5Momentum pin(parent.momentum()); double xs(partons[3]->mass()/pin.e()),xb(1.-xs); Lorentz5Momentum pspect(xs*pin),pdec(xb*pin); pspect.setMass(partons[3]->mass()); pdec.rescaleMass(); // Get the particle quark that is decaying long idQ, idSpec; idSpec = partons[3]->id(); idQ = (abs(parent.id())/1000)%10; if(!idQ) idQ = (abs(parent.id())/100)%10; // Now the odd case of a B_c where the c decays, not the b if(abs(idSpec) == idQ) idQ = (abs(parent.id())/10)%10; // change sign if spectator quark or antidiquark if((idSpec>0&&idSpec<6)||idSpec<-6) idQ = -idQ; // check if W products coloured bool Wcol = partons[0]->coloured(); // particle data object tcPDPtr dec = getParticleData(idQ); + // spin density matrix for the decaying quark + RhoDMatrix rhoin(PDT::Spin1Half,true); + if(parent.dataPtr()->iSpin()!=PDT::Spin0 && parent.spinInfo()) { + parent.spinInfo()->decay(); + RhoDMatrix rhoHadron = parent.spinInfo()->rhoMatrix(); + // particles with spin 0 diquark + if(abs(parent.id())==5122 || abs(parent.id())==4122 || + abs(parent.id())==5122 || abs(parent.id())==4122 || + abs(parent.id())==5132 || abs(parent.id())==4132 || + abs(parent.id())==5232 || abs(parent.id())==4232) { + for(unsigned int ix=0;ix<2;++ix) rhoin(ix,ix) = rhoHadron(ix,ix); + } + // particles with spin 1 diquark + else if(abs(parent.id())==5332 || abs(parent.id())==4332) { + rhoin(0,0) = 2./3.*rhoHadron(1,1)+1./3.*rhoHadron(0,0); + rhoin(1,1) = 2./3.*rhoHadron(0,0)+1./3.*rhoHadron(1,1); + } + } // momenta of the decay products vector pout(3,Lorentz5Momentum()); for(unsigned int ix=0;ix<3;++ix) pout[ix].setMass(partons[ix]->mass()); // charges of the exchanged boson and check if colour rearranged int c1 = dec ->iCharge()-partons[2]->dataPtr()->iCharge(); int c2 = partons[0]->dataPtr()->iCharge()+partons[1]->dataPtr()->iCharge(); bool rearranged = !(c1==c2&&abs(c1)==3); if(MECode==0) rearranged=false; if(rearranged) { int c3 = dec ->iCharge()-partons[1]->dataPtr()->iCharge(); int c4 = partons[0]->dataPtr()->iCharge()+partons[2]->dataPtr()->iCharge(); if(!(c3==c4&&abs(c3)==3)) { generator()->log() << "Unknown order for colour rearranged decay" << " in WeakPartonicDecayer::decay()\n"; generator()->log() << c1 << " " << c2 << " " << c3 << " " << c4 << "\n"; generator()->log() << parent << "\n" << dec->PDGName() << "\n"; for(unsigned int ix=0;ix<4;++ix) generator()->log() << *partons[ix] << "\n"; throw Exception() << "Unknown order for colour rearranged decay" << " in WeakPartonicDecayer::decay() " << Exception::runerror; } swap(pout[1] ,pout[2] ); swap(partons[1],partons[2]); } // decide if three or four body using prob bool threeBody = UseRandom::rnd() > _radprob; - // if four body not kinematically possible must be four body + // if four body not kinematically possible must be three body if(pdec.mass()constituentMass()+pout[0].mass()+ pout[1].mass()+pout[2].mass()) threeBody=true; // if code ==0 always three body if(MECode==0) threeBody=true; // three body decay if( threeBody ) { if(MECode==0) { Kinematics::threeBodyDecay(pdec,pout[1],pout[0],pout[2]); } else { // generate the kinematics double wgt(0.); Energy2 mb2max = sqr(pdec.mass() - pout[2].mass()); Energy2 mb2min = sqr(pout[0].mass() + pout[1].mass()); unsigned int ntry = 0; do { ++ntry; Energy2 mb2 = (mb2max-mb2min)*UseRandom::rnd()+mb2min; double CosAngle, AzmAngle; // perform first decay Lorentz5Momentum p01; p01.setMass(sqrt(mb2)); Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(pdec,pout[2].mass(),p01.mass(), CosAngle,AzmAngle,pout[2],p01); // perform second decay Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(p01,pout[0].mass(),pout[1].mass(), CosAngle,AzmAngle,pout[0],pout[1]); // kinematic piece of the weight wgt = Kinematics::pstarTwoBodyDecay(pdec.mass(),p01 .mass(),pout[2].mass())/pdec.mass()* Kinematics::pstarTwoBodyDecay(p01 .mass(),pout[0].mass(),pout[1].mass())/p01.mass(); - // piece to improve weight variation + // piece to improve weight variation (not kinematics dependent) wgt *= pdec.mass()/Kinematics::pstarTwoBodyDecay(pdec.mass(),sqrt(mb2min),pout[2].mass()); + // integration over m23^2 + wgt *= (mb2max-mb2min)/sqr(pdec.mass()); + // set momenta of particles + for(unsigned int ix=0;ixsetMomentum(pout[ix]); // matrix element piece - wgt *= 16.*(pdec*pout[1])*(pout[0]*pout[2])/sqr(mb2max-mb2min); + wgt *= threeBodyMatrixElement(dec,rhoin,pdec,partons); // check doesn't violate max if(wgt>_threemax) { ostringstream message; message << "Maximum weight for three-body decay " - << "violated in WeakPartonicDecayer2::decay()" + << "violated in WeakPartonicDecayer::decay()" << "Maximum = " << _threemax << " weight = " << wgt; generator()->logWarning( Exception(message.str(),Exception::warning) ); } } while( wgt < _threemax*UseRandom::rnd() && ntry < _maxtry ); if(ntry==_maxtry) throw Exception() << "Too many attempts to generate three body kinematics in " - << "WeakPartonicDecayer2::decay()" << Exception::eventerror; + << "WeakPartonicDecayer::decay()" << Exception::eventerror; } - // set momenta of particles - for(unsigned int ix=0;ixsetMomentum(pout[ix]); partons[3]->setMomentum(pspect); - // special for tau leptons to get correlations - if(abs(partons[0]->id())==ParticleID::tauminus|| - abs(partons[1]->id())==ParticleID::tauminus) - threeBodyMatrixElement(dec,pdec,partons); // set up the colour connections if(rearranged) swap(partons[1],partons[2]); if(Wcol) { if(partons[0]->data().iColour()==PDT::Colour3) partons[0]->antiColourNeighbour(partons[1]); else partons[0]-> colourNeighbour(partons[1]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[3]); } else { partons[2]-> colourNeighbour(partons[3]); } } // four body decay else { // generate the extra gluon partons.push_back(gluon->produceParticle()); partons.back()->set5Momentum(gluon->constituentMass()); // momentum of gluon pout.push_back(Lorentz5Momentum(gluon->constituentMass())); // generate the kinematics Energy2 ms2min(sqr(pout[0].mass()+pout[1].mass()+pout[2].mass())); Energy2 ms2max(sqr(pdec.mass()-pout[3].mass())); double wgt(0.); unsigned int ntry=0; bool initial = true; do { ++ntry; Energy2 ms2 = ms2min+UseRandom::rnd()*(ms2max-ms2min); Energy ms = sqrt(ms2); // and the W Energy2 mb2max = sqr(ms -pout[2].mass()); Energy2 mb2min = sqr(pout[0].mass()+pout[1].mass()); Energy2 mb2 = (mb2max-mb2min)*UseRandom::rnd()+mb2min; wgt = (mb2max-mb2min)/(ms2max-mb2min); // perform first decay Lorentz5Momentum ps; double CosAngle,AzmAngle; Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(pdec,pout[3].mass(),ms,CosAngle,AzmAngle,pout[3],ps); // generate the kinematics // perform second decay Kinematics::generateAngles(CosAngle,AzmAngle); Lorentz5Momentum p01; p01.setMass(sqrt(mb2)); Kinematics::twoBodyDecay(ps,pout[2].mass(),p01.mass(), CosAngle,AzmAngle,pout[2],p01); // perform third decay Kinematics::generateAngles(CosAngle,AzmAngle); Kinematics::twoBodyDecay(p01,pout[0].mass(),pout[1].mass(), CosAngle,AzmAngle,pout[0],pout[1]); // kinematic piece of the weight wgt *= 16.* Kinematics::pstarTwoBodyDecay(pdec.mass(),pout[3].mass(),ms )/pdec.mass()* Kinematics::pstarTwoBodyDecay(ms ,p01 .mass(),pout[2].mass())/ms* Kinematics::pstarTwoBodyDecay(p01 .mass(),pout[0].mass(),pout[1].mass())/p01.mass(); wgt *= fourBodyMatrixElement(pdec,pout[2],pout[0],pout[1],pout[3],Wcol,initial); // check doesn't violate max if(wgt>_threemax) { ostringstream message; message << "Maximum weight for four-body decay " << "violated in WeakPartonicDecayer::decay()" << "Maximum = " << _fourmax << " weight = " << wgt; generator()->logWarning( Exception(message.str(),Exception::warning) ); } } while ( wgt < _fourmax*UseRandom::rnd() && ntry < _maxtry ); if(ntry==_maxtry) throw Exception() << "Too many attempts to generate four body kinematics in " << "WeakPartonicDecayer::decay()" << Exception::eventerror; // set momenta of particles for(unsigned int ix=0;ix<3;++ix) partons[ix]->setMomentum(pout[ix]); partons[3]->setMomentum(pspect); partons[4]->setMomentum(pout[3]); // special for tau leptons to get correlations - if(abs(partons[0]->id())==ParticleID::tauminus|| - abs(partons[1]->id())==ParticleID::tauminus) - threeBodyMatrixElement(dec,pdec,partons); + threeBodyMatrixElement(dec,rhoin,pdec,partons); // set up the colour connections if(rearranged) swap(partons[1],partons[2]); // radiation from initial-state if(initial) { if(Wcol) { if(partons[0]->data().iColour()==PDT::Colour3) partons[0]->antiColourNeighbour(partons[1]); else partons[0]-> colourNeighbour(partons[1]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[4]); partons[3]-> colourNeighbour(partons[4]); } else { partons[2]-> colourNeighbour(partons[4]); partons[3]->antiColourNeighbour(partons[4]); } } // radiation from final-state else { if(partons[0]->data().iColour()==PDT::Colour3) { partons[0]->antiColourNeighbour(partons[4]); partons[1]-> colourNeighbour(partons[4]); } else { partons[0]-> colourNeighbour(partons[4]); partons[1]->antiColourNeighbour(partons[4]); } if(partons[2]->data().iColour()==PDT::Colour3) { partons[2]->antiColourNeighbour(partons[3]); } else { partons[2]-> colourNeighbour(partons[3]); } } } } return partons; } void WeakPartonicDecayer::persistentOutput(PersistentOStream & os) const { os << MECode << _radprob << _maxtry << _threemax << _fourmax; } void WeakPartonicDecayer::persistentInput(PersistentIStream & is, int) { is >> MECode >> _radprob >> _maxtry >> _threemax >> _fourmax; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigWeakPartonicDecayer("Herwig::WeakPartonicDecayer", "HwPartonicDecay.so"); void WeakPartonicDecayer::Init() { static ClassDocumentation documentation ("The WeakPartonicDecayer class performs partonic decays of hadrons containing a " "heavy quark."); static Switch interfaceMECode ("MECode", "The code for the type of matrix element to be used.", &WeakPartonicDecayer::MECode, 0, false, false); static SwitchOption interfaceMECodePhaseSpace (interfaceMECode, "PhaseSpace", "Phase space decays", 0); static SwitchOption interfaceMECodeWeak (interfaceMECode, "Weak", "Weak matrix element", 100); static Parameter interfaceRadiationProbability ("RadiationProbability", "The probability that QCD radiation produces an extra q qbar pair", &WeakPartonicDecayer::_radprob, 0., 0.0, 1., false, false, Interface::limited); static Parameter interfaceMaxTry ("MaxTry", "The maximum number of attempts to generate the kinematics", &WeakPartonicDecayer::_maxtry, 300, 10, 1000, false, false, Interface::limited); static Parameter interfaceThreeMax ("ThreeMax", "Maximum weight for sampling of three-body decays", &WeakPartonicDecayer::_threemax, 3.0, 1.0, 1000.0, false, false, Interface::limited); static Parameter interfaceFourMax ("FourMax", "Maximum weight for sampling of four-body decays", &WeakPartonicDecayer::_fourmax, 3.0, 1.0, 1000.0, false, false, Interface::limited); } double WeakPartonicDecayer::VAWt(Energy2 t0, Energy2 t1, Energy2 t2, InvEnergy4 t3) { return (t1-t0)*(t0-t2)*t3; } void WeakPartonicDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the PartonicDecayerBase base class PartonicDecayerBase::dataBaseOutput(output,false); output << "newdef " << name() << ":MECode " << MECode << " \n"; if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } -void WeakPartonicDecayer:: -threeBodyMatrixElement(tcPDPtr dec,Lorentz5Momentum & pdec, +double WeakPartonicDecayer:: +threeBodyMatrixElement(tcPDPtr dec, const RhoDMatrix & rhoin, + Lorentz5Momentum & pdec, ParticleVector & partons) const { // spinors LorentzSpinor w0[2],w2[2]; LorentzSpinorBar w1[2],w3[2]; // spinors for the decaying particle and first product if(dec->id()>0) { SpinorWaveFunction win(pdec,dec,0,incoming); w0[0] = win.dimensionedWave(); win.reset(1); w0[1] = win.dimensionedWave(); SpinorBarWaveFunction wout(partons[2]->momentum(), partons[2]->dataPtr(),0,outgoing); w1[0] = wout.dimensionedWave(); wout.reset(1); w1[1] = wout.dimensionedWave(); } else { SpinorBarWaveFunction win(pdec,dec,0,incoming); w1[0] = win.dimensionedWave(); win.reset(1); w1[1] = win.dimensionedWave(); SpinorWaveFunction wout(partons[2]->momentum(), partons[2]->dataPtr(),0,outgoing); w0[0] = wout.dimensionedWave(); wout.reset(1); w0[1] = wout.dimensionedWave(); } // spinors for the W decay products - bool taufirst = true; + bool lorder = true; if(partons[0]->id()<0) { SpinorWaveFunction wout2(partons[0]->momentum(), partons[0]->dataPtr(),0,outgoing); SpinorBarWaveFunction wout3(partons[1]->momentum(), partons[1]->dataPtr(),0,outgoing); + lorder = partons[0]->dataPtr()->charged(); w2[0] = wout2.dimensionedWave(); w3[0] = wout3.dimensionedWave(); wout2.reset(1); wout3.reset(1); w2[1] = wout2.dimensionedWave(); w3[1] = wout3.dimensionedWave(); - if(abs(partons[0]->id())!=ParticleID::tauminus) taufirst=false; } else { SpinorWaveFunction wout2(partons[1]->momentum(), partons[1]->dataPtr(),0,outgoing); SpinorBarWaveFunction wout3(partons[0]->momentum(), partons[0]->dataPtr(),0,outgoing); + lorder = partons[1]->dataPtr()->charged(); w2[0] = wout2.dimensionedWave(); w3[0] = wout3.dimensionedWave(); wout2.reset(1); wout3.reset(1); w2[1] = wout2.dimensionedWave(); w3[1] = wout3.dimensionedWave(); - if(abs(partons[1]->id())!=ParticleID::tauminus) taufirst=false; } + bool tau = abs(partons[0]->id())==ParticleID::tauminus || abs(partons[1]->id())==ParticleID::tauminus; // calculate the currents LorentzPolarizationVectorE Jbc[2][2],Jdec[2][2]; for(unsigned int ix=0;ix<2;++ix) { for(unsigned int iy=0;iy<2;++iy) { - Jbc [ix][iy] = w0[ix].leftCurrent(w1[iy]); - Jdec[ix][iy] = w2[ix].leftCurrent(w3[iy]); + if(dec->id()>0) + Jbc [ix][iy] = w0[ix].leftCurrent(w1[iy]); + else + Jbc [ix][iy] = w0[iy].leftCurrent(w1[ix]); + if(lorder) + Jdec[ix][iy] = w2[ix].leftCurrent(w3[iy]); + else + Jdec[ix][iy] = w2[iy].leftCurrent(w3[ix]); } } // compute the matrix element Complex me[2][2][2][2]; + double total=0.; for(unsigned int i0=0;i0<2;++i0) { for(unsigned int i1=0;i1<2;++i1) { for(unsigned int i2=0;i2<2;++i2) { for(unsigned int i3=0;i3<2;++i3) { me[i0][i1][i2][i3] = Jbc[i0][i1].dot(Jdec[i2][i3])/sqr(pdec.mass()); + total += rhoin(i0,i0).real()*norm(me[i0][i1][i2][i3]); } } } } - RhoDMatrix rho(PDT::Spin1Half); - for(unsigned int it1=0;it1<2;++it1) { - for(unsigned int it2=0;it2<2;++it2) { - for(unsigned int i0=0;i0<2;++i0) { - for(unsigned int i1=0;i1<2;++i1) { - for(unsigned int i2=0;i2<2;++i2) { - rho(it1,it2) += taufirst ? - me[i0][i1][it1][i2 ]*conj(me[i0][i1][it2][i2 ]) : - me[i0][i1][i2 ][it1]*conj(me[i0][i1][i2 ][it2]); + total *=2.; + if(tau) { + RhoDMatrix rho(PDT::Spin1Half); + for(unsigned int it1=0;it1<2;++it1) { + for(unsigned int it2=0;it2<2;++it2) { + for(unsigned int i0=0;i0<2;++i0) { + for(unsigned int i1=0;i1<2;++i1) { + for(unsigned int i2=0;i2<2;++i2) { + rho(it1,it2) += me[i0][i1][it1][i2 ]*conj(me[i0][i1][it2][i2 ]); + } } } } } + // normalize matrix to unit trace + rho.normalize(); + for(unsigned int ix=0;ix<2;++ix) { + if(abs(partons[ix]->id())!=ParticleID::tauminus) continue; + bool loc = partons[ix]->id() < 0; + // create the spin info object + FermionSpinPtr spin = new_ptr(FermionSpinInfo(partons[ix]->momentum(),true)); + // assign spinors + for(unsigned int iy=0;iy<2;++iy) { + spin->setBasisState(iy, loc ? w2[iy] : w3[iy].bar()); + } + // assign rho + spin->rhoMatrix() = rho; + // assign spin info + partons[ix]->spinInfo(spin); + } } - // normalize matrix to unit trace - rho.normalize(); - for(unsigned int ix=0;ix<2;++ix) { - if(abs(partons[ix]->id())!=ParticleID::tauminus) continue; - bool loc = partons[ix]->id() < 0; - // create the spin info object - FermionSpinPtr spin = new_ptr(FermionSpinInfo(partons[ix]->momentum(),true)); - // assign spinors - for(unsigned int iy=0;iy<2;++iy) { - spin->setBasisState(iy, loc ? w2[iy] : w3[iy].bar()); - } - // assign rho - spin->rhoMatrix() = rho; - // assign spin info - partons[ix]->spinInfo(spin); - } + return total; } double WeakPartonicDecayer:: fourBodyMatrixElement(Lorentz5Momentum & p0,Lorentz5Momentum & p1, Lorentz5Momentum & p2,Lorentz5Momentum & p3, Lorentz5Momentum & pg, bool Wcol, bool & initial) const { Energy2 d01(p0*p1),d02(p0*p2),d03(p0*p3),d0g(p0*pg); Energy2 d12(p1*p2),d13(p1*p3),d1g(p1*pg); Energy2 d23(p2*p3),d2g(p2*pg),d3g(p3*pg); Energy2 m02(sqr(p0.mass())),m12(sqr(p1.mass())),m22(sqr(p2.mass())), m32(sqr(p3.mass())); Energy2 mei = +1./d0g/d1g *( -d01*d12*d3g+d01*d03*d2g+2*d01*d03*d12 ) +1./d0g *( d12*d3g-d03*d12-d02*d03 ) +1./d1g *( d12*d13+d03*d2g+d03*d12 ) +m12/sqr(d1g)*( -d03*d2g-d03*d12 ) +m02/sqr(d0g)*( d12*d3g-d03*d12 ); Energy2 mef = !Wcol ? ZERO : +1./d2g/d3g *( d0g*d12*d23+d03*d1g*d23+2*d03*d12*d23 ) +1./d2g *( d03*d1g+d03*d12-d02*d12 ) +1./d3g *( d0g*d12-d03*d13+d03*d12 ) +m32/sqr(d3g)*( -d0g*d12-d03*d12 ) +m22/sqr(d2g)*( -d03*d1g-d03*d12 ); initial = mef/(mei+mef)