diff --git a/Decay/Baryon/Baryon1MesonDecayerBase.cc b/Decay/Baryon/Baryon1MesonDecayerBase.cc --- a/Decay/Baryon/Baryon1MesonDecayerBase.cc +++ b/Decay/Baryon/Baryon1MesonDecayerBase.cc @@ -1,976 +1,976 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the Baryon1MesonDecayerBase class. // #include "Baryon1MesonDecayerBase.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" using namespace Herwig; using namespace ThePEG::Helicity; AbstractNoPIOClassDescription Baryon1MesonDecayerBase::initBaryon1MesonDecayerBase; // Definition of the static class description member. 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."); } // 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 & inpart, const ParticleVector & decay, MEOption meopt) const { double me(0.); // decide which matrix element we are doing // incoming spin-1/2 particle if(inpart.dataPtr()->iSpin()==2) { // decay to spin-1/2 particle if(decay[0]->dataPtr()->iSpin()==2) { // scalar meson if(decay[1]->dataPtr()->iSpin()==1) me=halfHalfScalar(ichan,inpart,decay,meopt); // vector meson else if(decay[1]->dataPtr()->iSpin()==3) me=halfHalfVector(ichan,inpart,decay,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // decay to spin-3/2 particle else if(decay[0]->dataPtr()->iSpin()==4) { // scalar meson if(decay[1]->dataPtr()->iSpin()==1) me=halfThreeHalfScalar(ichan,inpart,decay,meopt); // vector meson else if(decay[1]->dataPtr()->iSpin()==3) me=halfThreeHalfVector(ichan,inpart,decay,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(inpart.dataPtr()->iSpin()==4) { // decay to spin-1/2 particle if(decay[0]->dataPtr()->iSpin()==2) { // scalar meson if(decay[1]->dataPtr()->iSpin()==1) me=threeHalfHalfScalar(ichan,inpart,decay,meopt); // vector meson else if(decay[1]->dataPtr()->iSpin()==3) me=threeHalfHalfVector(ichan,inpart,decay,meopt); else throw DecayIntegratorError() << "Unknown outgoing meson spin in " << "Baryon1MesonDecayerBase::me2()" << Exception::abortnow; } // decay to spin-3/2 particle else if(decay[0]->dataPtr()->iSpin()==4) { // scalar meson if(decay[1]->dataPtr()->iSpin()==1) me=threeHalfThreeHalfScalar(ichan,inpart,decay,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 & inpart, const ParticleVector & decay,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(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); return 0.; } // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,decay[0],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf,decay[0],outgoing); } // get the couplings Complex A,B; halfHalfScalarCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(),A,B); Complex left,right,meout; // coupling for an incoming particle if(inpart.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(decay[0]->id()>0){ispin[0]=iy;ispin[1]=ix;} else{ispin[0]=ix;ispin[1]=iy;} - (*ME())(ispin)=_inHalf[iy].generalScalar(_inHalfBar[ix],left,right)/inpart.mass(); + (*ME())(ispin)=Complex(_inHalf[iy].generalScalar(_inHalfBar[ix],left,right)/inpart.mass()); // output += norm(ME()(ispin)); } } // test of the matrix elemen// t // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[1]->mass()); // Energy Qp(sqrt(sqr(m1+m2)-sqr(m3))),Qm(sqrt(sqr(m1-m2)-sqr(m3))); // Complex h1(2.*Qp*A/inpart.mass()),h2(-2.*Qm*B/inpart.mass()); // generator()->log() << "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; // generator()->log() << "testing alpha " << // (norm(0.5*(h1+h2))-norm(0.5*(h1-h2)))/ // (norm(0.5*(h1+h2))+norm(0.5*(h1-h2))) << "\n"; // Energy pcm(Kinematics::pstarTwoBodyDecay(m1,m2,m3)); // generator()->log() << "testing masses " << m1/GeV << " " << m2/GeV << " " << m3/GeV // << "\n"; // generator()->log() << "testing partial " << pcm*0.5*output/8./Constants::pi/MeV // << "\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 & inpart, const ParticleVector & decay,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=decay[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } VectorWaveFunction::constructSpinInfo(_inVec,decay[1],outgoing,true,photon); return 0.; } // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,decay[0],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf,decay[0],outgoing); } VectorWaveFunction::calculateWaveFunctions(_inVec,decay[1],outgoing,photon); // get the couplings Complex A1,A2,B1,B2; halfHalfVectorCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(), A1,A2,B1,B2); Complex lS,rS,lV,rV; complex scalar; // couplings for an incoming particle if(inpart.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(inpart.mass()+decay[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(decay[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(inpart.momentum())/msum; (*ME())(ispin)=(svec.dot(_inVec[ispin[2]])+prod*scalar)/inpart.mass(); // output += norm(ME()(ispin)); } } } // test of the matrix element // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[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)); // Complex h1(2.*r2*Qp*B1/inpart.mass()),h2(-2.*r2*Qm*A1/inpart.mass()), // h3(2./m3*(Qp*(m1-m2)*B1-Qm*m1*B2*pcm/(m1+m2))/inpart.mass()), // h4(2./m3*(Qm*(m1+m2)*A1+Qp*m1*A2*pcm/(m1+m2))/inpart.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"; // 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 & inpart, const ParticleVector & decay, 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(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); return 0.; } // spinors for the decay product LorentzPolarizationVector in=UnitRemoval::InvE*inpart.momentum(); if(inpart.id()>0) { RSSpinorBarWaveFunction:: calculateWaveFunctions(_inThreeHalfBar,decay[0],outgoing); _inHalfBar.resize(_inThreeHalfBar.size()); for(unsigned int ix=0;ix<_inThreeHalfBar.size();++ix) _inHalfBar[ix] = _inThreeHalfBar[ix].dot(in); } else { RSSpinorWaveFunction:: calculateWaveFunctions(_inThreeHalf,decay[0],outgoing); _inHalf.resize(_inThreeHalf.size()); for(unsigned int ix=0;ix<_inThreeHalf.size();++ix) _inHalf[ix] = _inThreeHalf[ix].dot(in); } // get the couplings Complex A,B,left,right; Energy msum(inpart.mass()+decay[0]->mass()); halfThreeHalfScalarCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(), A,B); // incoming particle if(inpart.id()>0) { left=(A-B); right=(A+B); } // incoming anti-particle else { left=conj(A+B); right=conj(A-B); } vector ispin(3,0); //Complex output(0.); for(unsigned ixa=0;ixa<2;++ixa) { for(unsigned int iya=0;iya<4;++iya) { unsigned int ix(iya),iy(ixa); if(decay[0]->id()<0) swap(ix,iy); ispin[0]=ixa; ispin[1]=iya; complex value = _inHalf[iy].generalScalar(_inHalfBar[ix],left,right) *UnitRemoval::E/inpart.mass()/msum; (*ME())(ispin) = value; //output+= norm(ME()(ispin)); } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[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 " << inpart.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2))/sqr(inpart.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2))/sqr(inpart.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 & inpart, const ParticleVector & decay, 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=decay[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(inpart.id()>0) SpinorWaveFunction ::calculateWaveFunctions(_inHalf,_rho, const_ptr_cast(&inpart), incoming); else SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,_rho, const_ptr_cast(&inpart), incoming); // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { SpinorWaveFunction:: constructSpinInfo(_inHalf,const_ptr_cast(&inpart),incoming,true); RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else { SpinorBarWaveFunction:: constructSpinInfo(_inHalfBar,const_ptr_cast(&inpart),incoming,true); RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } VectorWaveFunction::constructSpinInfo(_inVec,decay[1],outgoing,true,photon); return 0.; } LorentzPolarizationVector in=UnitRemoval::InvE*inpart.momentum(); if(inpart.id()>0) { RSSpinorBarWaveFunction:: calculateWaveFunctions(_inThreeHalfBar,decay[0],outgoing); _inHalfBar.resize(_inThreeHalfBar.size()); for(unsigned int ix=0;ix<_inThreeHalfBar.size();++ix) _inHalfBar[ix] = _inThreeHalfBar[ix].dot(in); } else { RSSpinorWaveFunction:: calculateWaveFunctions(_inThreeHalf,decay[0],outgoing); _inHalf.resize(_inThreeHalf.size()); for(unsigned int ix=0;ix<_inThreeHalf.size();++ix) _inHalf[ix] = _inThreeHalf[ix].dot(in); } ME()->zero(); VectorWaveFunction::calculateWaveFunctions(_inVec,decay[1],outgoing,photon); // get the couplings Complex A1,A2,A3,B1,B2,B3; halfThreeHalfVectorCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(), A1,A2,A3,B1,B2,B3); Energy msum(inpart.mass()+decay[0]->mass()); Complex lS,rS,lV,rV,left,right; // incoming particle if(inpart.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(decay[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(inpart.momentum())/msum; (*ME())(ispin) += (svec.dot(_inVec[iz])+prod*scalar)* UnitRemoval::E/msum/inpart.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(decay[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(decay[0]->id()>0) stemp = _inHalf[ixa]; else sbtemp = _inHalfBar[ixa]; - (*ME())(ispin) += stemp.generalScalar(sbtemp,left,right)/inpart.mass(); + (*ME())(ispin) += Complex(stemp.generalScalar(sbtemp,left,right)/inpart.mass()); } } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[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 " << inpart.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(inpart.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(inpart.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 & inpart, const ParticleVector & decay, 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(inpart.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&inpart), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&inpart), incoming); } // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { RSSpinorWaveFunction:: constructSpinInfo(_inThreeHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { RSSpinorBarWaveFunction:: constructSpinInfo(_inThreeHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); return 0.; } // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,decay[0],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf,decay[0],outgoing); } LorentzPolarizationVector out=UnitRemoval::InvE*decay[0]->momentum(); if(inpart.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=inpart.mass()+decay[0]->mass(); threeHalfHalfScalarCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(), A,B); Complex left,right; // incoming particle if(inpart.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(decay[0]->id()<0) swap(ix,iy); ispin[0]=iya; ispin[1]=ixa; - (*ME())(ispin) = _inHalf[iy].generalScalar(_inHalfBar[ix],left,right)* - UnitRemoval::E/msum/inpart.mass(); + (*ME())(ispin) = Complex(_inHalf[iy].generalScalar(_inHalfBar[ix],left,right)* + UnitRemoval::E/msum/inpart.mass()); } } double output = (ME()->contract(_rho)).real(); // test of the matrix element // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[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 " << inpart.id() << " " // << output << " " // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(inpart.mass()) << " " // << 0.125*(h1*conj(h1)+h2*conj(h2))/sqr(inpart.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 & inpart, const ParticleVector & decay, 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(inpart.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&inpart), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&inpart), incoming); } // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { RSSpinorWaveFunction:: constructSpinInfo(_inThreeHalf,const_ptr_cast(&inpart),incoming,true); RSSpinorBarWaveFunction::constructSpinInfo(_inThreeHalfBar, decay[0],outgoing,true); } else { RSSpinorBarWaveFunction:: constructSpinInfo(_inThreeHalfBar,const_ptr_cast(&inpart),incoming,true); RSSpinorWaveFunction::constructSpinInfo(_inThreeHalf, decay[0],outgoing,true); } ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true); return 0.; } // spinors for the decay product LorentzPolarizationVector in=UnitRemoval::InvE*inpart.momentum(); if(inpart.id()>0) { RSSpinorBarWaveFunction:: calculateWaveFunctions(_inThreeHalfBar,decay[0],outgoing); } else { RSSpinorWaveFunction:: calculateWaveFunctions(_inThreeHalf,decay[0],outgoing); } _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(inpart.mass()+decay[0]->mass()); threeHalfThreeHalfScalarCoupling(imode(),inpart.mass(),decay[0]->mass(), decay[1]->mass(),A1,A2,B1,B2); Complex left1,right1,left2,right2; // incoming particle if(inpart.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(decay[0]->id()<0) swap(ix,iy); ispin[0]=iya; ispin[1]=ixa; - (*ME())(ispin)=(_inThreeHalf[iy].generalScalar(_inThreeHalfBar[ix],left1,right1) - +_inHalf[iy].generalScalar( _inHalfBar[ix],left2,right2) - *UnitRemoval::E2/sqr(msum))/inpart.mass(); + (*ME())(ispin)=Complex((_inThreeHalf[iy].generalScalar(_inThreeHalfBar[ix],left1,right1) + +_inHalf[iy].generalScalar( _inHalfBar[ix],left2,right2) + *UnitRemoval::E2/sqr(msum))/inpart.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 & inpart, const ParticleVector & decay,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=decay[1]->id()==ParticleID::gamma; // spinors etc for the decaying particle if(meopt==Initialize) { // spinors and rho if(inpart.id()>0) { RSSpinorWaveFunction ::calculateWaveFunctions(_inThreeHalf,_rho, const_ptr_cast(&inpart), incoming); } else { RSSpinorBarWaveFunction::calculateWaveFunctions(_inThreeHalfBar,_rho, const_ptr_cast(&inpart), incoming); } // matrix element } // setup spin info when needed if(meopt==Terminate) { // for the decaying particle if(inpart.id()>0) { RSSpinorWaveFunction:: constructSpinInfo(_inThreeHalf,const_ptr_cast(&inpart),incoming,true); SpinorBarWaveFunction::constructSpinInfo(_inHalfBar,decay[0],outgoing,true); } else { RSSpinorBarWaveFunction:: constructSpinInfo(_inThreeHalfBar,const_ptr_cast(&inpart),incoming,true); SpinorWaveFunction::constructSpinInfo(_inHalf,decay[0],outgoing,true); } VectorWaveFunction::constructSpinInfo(_inVec,decay[1],outgoing,true,photon); return 0.; } // spinors for the decay product if(inpart.id()>0) { SpinorBarWaveFunction::calculateWaveFunctions(_inHalfBar,decay[0],outgoing); } else { SpinorWaveFunction ::calculateWaveFunctions(_inHalf,decay[0],outgoing); } LorentzPolarizationVector out=UnitRemoval::InvE*decay[0]->momentum(); if(inpart.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); } ME()->zero(); VectorWaveFunction::calculateWaveFunctions(_inVec,decay[1],outgoing,photon); // get the couplings Complex A1,A2,A3,B1,B2,B3,prod,meout; threeHalfHalfVectorCoupling(imode(),inpart.mass(),decay[0]->mass(),decay[1]->mass(), A1,A2,A3,B1,B2,B3); Energy msum(inpart.mass()+decay[0]->mass()); Complex lS,rS,lV,rV,left,right; // incoming particle if(inpart.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(decay[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(decay[0]->momentum())/msum; (*ME())(ispin) += (svec.dot(_inVec[iz])+prod*scalar)* UnitRemoval::E/msum/inpart.mass(); } } // the piece where the vector spinor is dotted with the polarization vector for(unsigned iz=0;iz<3;++iz) { ispin[2]=iz; if(decay[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(decay[0]->id()>0) sbtemp = _inHalfBar[ixa]; else stemp = _inHalf[ixa]; - (*ME())(ispin) += stemp.generalScalar(sbtemp,left,right)/inpart.mass(); + (*ME())(ispin) += Complex(stemp.generalScalar(sbtemp,left,right)/inpart.mass()); } } } double output = (ME()->contract(_rho)).real(); // testing code // Energy m1(inpart.mass()),m2(decay[0]->mass()),m3(decay[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(-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)); // complex 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 " << inpart.id() << " " // << output << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(inpart.mass()) << " " // << 0.25*(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+ // h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/sqr(inpart.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/MatrixElement/Matchbox/Utility/AmplitudeCache.h b/MatrixElement/Matchbox/Utility/AmplitudeCache.h --- a/MatrixElement/Matchbox/Utility/AmplitudeCache.h +++ b/MatrixElement/Matchbox/Utility/AmplitudeCache.h @@ -1,346 +1,346 @@ // -*- C++ -*- // // AmplitudeCache.h 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. // #ifndef HERWIG_AmplitudeCache_H #define HERWIG_AmplitudeCache_H #include "Herwig/MatrixElement/Matchbox/Utility/SpinorHelicity.h" #include "ThePEG/Config/algorithm.h" #include namespace Herwig { using namespace ThePEG; using std::array; namespace SpinorHelicity { /** * \ingroup Matchbox * \author Simon Platzer * * \brief Caching for amplitudes using spinor helicity techniques. * */ template class AmplitudeCache { typedef map > AmplitudeCacheMap; typedef map > > CurrentCacheMap; /** * Maximum N we can handle, SYM_N is storage size for a symmetric matrix of N * N elements */ enum { MAX_N = 7, SYM_N = MAX_N*(MAX_N+1)/2 }; /** * The number of points */ int theNPoints; /** * The energy scale to obtain dimensionless * quantities. */ mutable Energy theScale; /** * Masses indexed by partons */ mutable array theMasses; /** * Momenta indexed by partons */ mutable array theMomenta; /** * Crossing signs indexed by partons */ mutable array theCrossingSigns; /** * Plus spinors indexed by partons */ mutable array thePlusSpinors; /** * Plus conjugate spinors indexed by partons */ mutable array thePlusConjugateSpinors; /** * Invariants indexed by partons */ mutable array theInvariants; /** * Flag products to be recalculated */ mutable array getInvariant; /** * Spinor products indexed by partons */ mutable array thePlusProducts; /** * Flag products to be recalculated */ mutable array getPlusProduct; /** * Spinor currents indexed by partons */ mutable array,SYM_N> thePlusCurrents; /** * Flag currents to be recalculated */ mutable array getPlusCurrent; /** * Cache intermediate amplitudes by index */ mutable AmplitudeCacheMap theCachedAmplitudes; /** * The last query for a cached amplitude */ mutable typename AmplitudeCacheMap::iterator theLastAmplitude; /** * Cache intermediate currents by index */ mutable CurrentCacheMap theCachedCurrents; /** * The last query for a cached current */ mutable typename CurrentCacheMap::iterator theLastCurrent; /** * Helper function to index symmetric arrays, assumes i <= j. * Usual indexing function (N*i + j) corrected by triangle number for i-th row. */ inline size_t idx(size_t i, size_t j) const { return MAX_N * i - (i + 1) * i / 2 + j; } /** * Helper to reset flags */ struct boolResetter { void operator()(pair >& flag) const { flag.second.first = true; } void operator()(pair > >& flag) const { flag.second.first = true; } }; public: /** * Constructor */ AmplitudeCache() : theNPoints(0) {} /** * Prepare for n-point amplitude */ void nPoints(int n); /** * Return the number of points */ int nPoints() const { return theNPoints; } /** * Set the energy scale to obtain dimensionless * quantities and flag all quantities to be recalculated. */ void amplitudeScale(Energy s) const; /** * Set the momentum for the k'th parton * and its associated mass. */ void momentum(int k, const LorentzMomentum& p, bool getSpinors = true, Energy mass = ZERO) const; /** * Reset flags */ void reset() const; public: /** * Return the momentum for the k'th parton */ LorentzVector momentum(int k) const { return theMomenta[k]/theScale; } /** * Get the energy scale to obtain dimensionless * quantities and flag all quantities to be recalculated. */ Energy amplitudeScale() const { return theScale; } /** * Return the mass associated to the k'th parton */ double mass(int k) const { return theMasses[k]; } /** * Return the crossing sign for the * i'th parton */ int crossingSign(int i) const { return theCrossingSigns[i]; } /** * Return the crossing sign for the * i'th and j'th parton */ double crossingSign(int i, int j) const { return theCrossingSigns[i]*theCrossingSigns[j]; } /** * Return (ij) */ double invariant(int i, int j) const { if ( i == j ) return 0.; if ( i > j ) swap(i,j); if ( getInvariant[idx(i,j)] ) { getInvariant[idx(i,j)] = false; theInvariants[idx(i,j)] = 2.*(momentum(i)*momentum(j)); } return theInvariants[idx(i,j)]; } /** * Return */ Complex plusProduct(int i, int j) const { if ( i== j ) return 0.; bool swapij = (i > j); if ( swapij ) swap(i,j); if ( getPlusProduct[idx(i,j)] ) { getPlusProduct[idx(i,j)] = false; thePlusProducts[idx(i,j)] = - PlusSpinorProduct(thePlusConjugateSpinors[i], - thePlusSpinors[j]).eval() / theScale; + Complex(PlusSpinorProduct(thePlusConjugateSpinors[i], + thePlusSpinors[j]).eval() / theScale); } return swapij ? -thePlusProducts[idx(i,j)] : thePlusProducts[idx(i,j)]; } /** * Return [ij] */ Complex minusProduct(int i, int j) const { if ( i== j ) return 0.; return -crossingSign(i,j)*conj(plusProduct(i,j)); } /** * Return plusCurrent(int i, int j) const { bool swapij = (i > j); if ( swapij ) swap(i,j); if ( getPlusCurrent[idx(i,j)] ) { getPlusCurrent[idx(i,j)] = false; if ( i != j ) { thePlusCurrents[idx(i,j)] = PlusSpinorCurrent(thePlusConjugateSpinors[i], MinusSpinor(theMomenta[j])).eval() / theScale; } else { thePlusCurrents[idx(i,j)] = 2.*momentum(i); } } return swapij ? crossingSign(i,j)*thePlusCurrents[idx(i,j)].conjugate() : thePlusCurrents[idx(i,j)]; } /** * Return [i|\gamma^\mu|j> */ LorentzVector minusCurrent(int i, int j) const { return plusCurrent(j,i); } public: /** * Return true, if the given amplitude * needs to be recalculated. */ bool getAmplitude(const AmplitudeKey& key) const { static Complex czero; if ( ( theLastAmplitude = theCachedAmplitudes.find(key) ) == theCachedAmplitudes.end() ) { theLastAmplitude = theCachedAmplitudes.insert(make_pair(key,make_pair(true,czero))).first; } return theLastAmplitude->second.first; } /** * Cache an amplitude */ void cacheAmplitude(Complex amp) const { theLastAmplitude->second = make_pair(false,amp); } /** * Return a cached amplitude */ const Complex& cachedAmplitude() const { return theLastAmplitude->second.second; } /** * Return true, if the given current * needs to be recalculated. */ bool getCurrent(const AmplitudeKey& key) const { static LorentzVector czero; if ( ( theLastCurrent = theCachedCurrents.find(key) ) == theCachedCurrents.end() ) { theLastCurrent = theCachedCurrents.insert(make_pair(key,make_pair(true,czero))).first; } return theLastCurrent->second.first; } /** * Cache an current */ void cacheCurrent(const LorentzVector& curr) const { theLastCurrent->second = make_pair(false,curr); } /** * Return a cached current */ const LorentzVector& cachedCurrent() const { return theLastCurrent->second.second; } }; } } #include "AmplitudeCache.tcc" #endif // HERWIG_AmplitudeCache_H diff --git a/Models/StandardModel/SMHPPVertex.cc b/Models/StandardModel/SMHPPVertex.cc --- a/Models/StandardModel/SMHPPVertex.cc +++ b/Models/StandardModel/SMHPPVertex.cc @@ -1,282 +1,282 @@ // -*- C++ -*- // // SMHPPVertex.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 SMHPPVertex class. // #include "SMHPPVertex.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/Looptools/clooptools.h" using namespace Herwig; using namespace ThePEG; void SMHPPVertex::persistentOutput(PersistentOStream & os) const { os << _theSM << ounit(_mw,GeV) << massopt << _minloop << _maxloop << _CoefRepresentation; } void SMHPPVertex::persistentInput(PersistentIStream & is, int) { is >> _theSM >> iunit(_mw, GeV) >> massopt >> _minloop >> _maxloop >> _CoefRepresentation; } ClassDescription SMHPPVertex::initSMHPPVertex; // Definition of the static class description member. void SMHPPVertex::Init() { static ClassDocumentation documentation ("This class implements the h0->gamma,gamma vertex."); static Parameter interfaceMinQuarkInLoop ("MinQuarkInLoop", "The minimum flavour of the quarks to include in the loops", &SMHPPVertex::_minloop, 6, 1, 6, false, false, Interface::limited); static Parameter interfaceMaxQuarkInLoop ("MaxQuarkInLoop", "The maximum flavour of the quarks to include in the loops", &SMHPPVertex::_maxloop, 6, 1, 6, false, false, Interface::limited); static Switch interfaceMassOption ("LoopMassScheme", "Switch for the treatment of the masses in the loops ", &SMHPPVertex::massopt, 2, false, false); static SwitchOption interfaceHeavyMass (interfaceMassOption, "PoleMasses", "The loop is calculcated with the pole quark masses", 1); static SwitchOption interfaceNormalMass (interfaceMassOption, "RunningMasses", "running quark masses are taken in the loop", 2); static Switch interfaceScheme ("CoefficientScheme", "Which scheme for the tensor coefficients is applied", &SMHPPVertex::_CoefRepresentation, 1, false, false); static SwitchOption interfaceSchemeSimplified (interfaceScheme, "Simplified", "Represection suitable for the simplified the H-g-g and H-gamma-gamma vertices", 1); static SwitchOption interfaceSchemeGeneral (interfaceScheme, "General", "Represection suitable for the Passarino-Veltman tensor reduction scheme", 2); } void SMHPPVertex::setCoupling(Energy2 q2, tcPDPtr part2, tcPDPtr part3, tcPDPtr part1) { assert( part1->id() == ParticleID::h0 && part2->id() == ParticleID::gamma && part3->id() == ParticleID::gamma ); int Qminloop = _minloop; int Qmaxloop = _maxloop; if (_maxloop < _minloop) { Qmaxloop=_minloop; Qminloop=_maxloop; } switch (_CoefRepresentation) { case 1: { if(q2 != _q2last||_couplast==0.) { double g = weakCoupling(q2); double e2 = sqr(electroMagneticCoupling(q2)); _couplast = UnitRemoval::E * e2 * g / 8. / _mw/ sqr(Constants::pi); _q2last = q2; } norm(_couplast); Complex loop(0.); // quark loops for ( int i = Qminloop; i <= Qmaxloop; ++i ) { tcPDPtr qrk = getParticleData(i); Energy mass = (2 == massopt) ? _theSM->mass(q2,qrk) : qrk->mass(); Charge charge = qrk->charge(); - loop += 3.*sqr(charge/ThePEG::Units::eplus) * Af(sqr(mass)/q2); + loop += Complex(3.*sqr(charge/ThePEG::Units::eplus) * Af(sqr(mass)/q2)); } // lepton loops int Lminloop = 3; // still fixed value int Lmaxloop = 3; // still fixed value for (int i = Lminloop; i <= Lmaxloop; ++i) { tcPDPtr lpt = getParticleData(9 + 2*i); Energy mass = (2 == massopt) ? _theSM->mass(q2,lpt) : lpt->mass(); Charge charge = lpt->charge(); - loop += sqr(charge/ThePEG::Units::eplus) * Af(sqr(mass)/q2); + loop += Complex(sqr(charge/ThePEG::Units::eplus) * Af(sqr(mass)/q2)); } // W loop loop += Aw(sqr(_mw)/q2); a00(loop); a11(0.0); a12(0.0); a21(-loop); a22(0.0); aEp(0.0); break; } case 2: { if(q2 != _q2last||_couplast==0.) { Looptools::clearcache(); double e = electroMagneticCoupling(q2); _couplast = pow(e,3)/sqrt(sin2ThetaW()); _q2last = q2; } norm(_couplast); // quarks int delta = Qmaxloop - Qminloop + 1; type.resize(delta,PDT::SpinUnknown); masses.resize(delta,ZERO); for (int i = 0; i < delta; ++i) { tcPDPtr q = getParticleData(_minloop+i); type[i] = PDT::Spin1Half; masses[i] = (2 == massopt) ? _theSM->mass(q2,q) : q->mass(); double copl = -masses[i]*3.*sqr(q->iCharge()/3.)/_mw/2.; couplings.push_back(make_pair(copl, copl)); } // tau type.push_back(PDT::Spin1Half); tcPDPtr tau = getParticleData(ParticleID::tauminus); masses.push_back(_theSM->mass(q2,tau)); double copl = -masses.back()*sqr(tau->iCharge()/3.)/_mw/2.; couplings.push_back(make_pair(copl, copl)); // W type.push_back(PDT::Spin1); masses.push_back(_mw); const double mw = UnitRemoval::InvE*_mw; couplings.push_back(make_pair(mw,mw)); setNParticles(delta+2); VVSLoopVertex::setCoupling(q2, part1, part2, part3); break; } } } Complex SMHPPVertex::Af(const double tau) const { return tau*(4. - W2(tau)*(1. - 4.*tau)); } Complex SMHPPVertex::Aw(const double tau) const { return 0.5*(-3.*W2(tau)*tau*(4.*tau - 2.) - 12.*tau - 2.); } Complex SMHPPVertex::W2(double lambda) const { double pi = Constants::pi; if (0.0 == lambda) return 0.0; if (lambda < 0.0) return 4.*sqr(asinh(0.5*sqrt(-1./lambda))); double root(0.5*sqrt(1./lambda)); Complex ac(0.); // formulae from NPB297,221 if(root < 1.) { ac = -sqr(asin(root)); } else { double ex = acosh(root); ac = sqr(ex) - 0.25*sqr(pi) - pi*ex*Complex(0.,1.); } return 4.*ac; } SMHPPVertex::SMHPPVertex() :_couplast(0.),_q2last(),_mw(),massopt(1), _minloop(6),_maxloop(6),_CoefRepresentation(1) { orderInGs(0); orderInGem(3); } // functions for loops for testing // namespace { // Complex F0(double tau) { // Complex ft; // if(tau>=1.) // ft = sqr(asin(1./sqrt(tau))); // else { // double etap = 1.+sqrt(1.-tau); // double etam = 1.-sqrt(1.-tau); // ft = -0.25*sqr(log(etap/etam)-Constants::pi*Complex(0.,1.)); // } // return tau*(1.-tau*ft); // } // Complex FHalf(double tau,double eta) { // Complex ft; // if(tau>=1.) // ft = sqr(asin(1./sqrt(tau))); // else { // double etap = 1.+sqrt(1.-tau); // double etam = 1.-sqrt(1.-tau); // ft = -0.25*sqr(log(etap/etam)-Constants::pi*Complex(0.,1.)); // } // return -2.*tau*(eta+(1.-tau*eta)*ft); // } // Complex F1(double tau) { // Complex ft; // if(tau>=1.) // ft = sqr(asin(1./sqrt(tau))); // else { // double etap = 1.+sqrt(1.-tau); // double etam = 1.-sqrt(1.-tau); // ft = -0.25*sqr(log(etap/etam)-Constants::pi*Complex(0.,1.)); // } // return 2.+3.*tau+3.*tau*(2.-tau)*ft; // } // } void SMHPPVertex::doinit() { //PDG codes for particles at vertices addToList(22,22,25); _theSM = dynamic_ptr_cast(generator()->standardModel()); if( !_theSM ) throw InitException() << "SMHGGVertex::doinit() - The pointer to the SM object is null." << Exception::abortnow; _mw = getParticleData(ThePEG::ParticleID::Wplus)->mass(); VVSLoopVertex::doinit(); // // code to test the partial width // Energy mh = getParticleData(25)->mass(); // Complex I(0.); // for(long ix=int(_minloop);ix<=int(_maxloop);++ix) { // tcPDPtr qrk = getParticleData(ix); // Energy mt = (2 == massopt) ? _theSM->mass(sqr(mh),qrk) : qrk->mass(); // double tau = sqr(2.*mt/mh); // I += 3.*sqr(double(qrk->iCharge())/3.)*FHalf(tau,1.); // cerr << "testing half " << FHalf(tau,1) << " " << Af(0.25*tau) << "\n"; // } // for(long ix=15;ix<=15;++ix) { // tcPDPtr qrk = getParticleData(ix); // Energy mt = (2 == massopt) ? _theSM->mass(sqr(mh),qrk) : qrk->mass(); // double tau = sqr(2.*mt/mh); // I += sqr(double(qrk->iCharge())/3.)*FHalf(tau,1.); // } // I += F1(sqr(2.*_mw/mh)); // Energy width = sqr(weakCoupling(sqr(mh))*sqr(electroMagneticCoupling(sqr(mh)))) // /1024./pow(Constants::pi,5)/16.*sqr(mh/_mw)*mh*std::norm(I); // cerr << "testing anal " << width/GeV << "\n"; if(loopToolsInitialized()) Looptools::ltexi(); }