diff --git a/Decay/Baryon/Baryon1MesonDecayerBase.cc b/Decay/Baryon/Baryon1MesonDecayerBase.cc --- a/Decay/Baryon/Baryon1MesonDecayerBase.cc +++ b/Decay/Baryon/Baryon1MesonDecayerBase.cc @@ -1,978 +1,978 @@ // -*- 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" 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."); } // 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/Decay/FormFactors/ISGWFormFactor.h b/Decay/FormFactors/ISGWFormFactor.h --- a/Decay/FormFactors/ISGWFormFactor.h +++ b/Decay/FormFactors/ISGWFormFactor.h @@ -1,279 +1,280 @@ // -*- C++ -*- // // ISGWFormFactor.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_ISGWFormFactor_H #define HERWIG_ISGWFormFactor_H // // This is the declaration of the & h, Complex & k, complex & bp, complex & bm) const; //@} /** * Output the setup information for the particle database * @param os The stream to output the information to * @param header Whether or not to output the information for MySQL * @param create Whether or not to add a statement creating the object */ virtual void dataBaseOutput(ofstream & os,bool header,bool create) const; public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** The member which implements all the different form-factors * @param q2 The scale \f$q^2\f$. * @param iloc The location in the form-factor list. * @param id0 The PDG code of the incoming meson. * @param id1 The PDG code of the outgoing meson. * @param m0 The mass of the incoming meson. * @param m1 The mass of the outgoing meson. * @param f1 The first form-factor. * @param f2 The second form-factor. * @param f3 The third form-factor. * @param f4 The fourth form-factor. */ void formFactor(Energy2 q2,unsigned int iloc,int id0,int id1,Energy m0, Energy m1,Complex & f1,Complex & f2, Complex & f3,Complex & f4) const; // general member to calculate all the form-factors protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** + * Private and non-existent assignment operator. */ ISGWFormFactor & operator=(const ISGWFormFactor &) = delete; private: /** * The relativistic compensation factor */ double _kappa; /** @name Quark masses */ //@{ /** * The down quark mass */ Energy _mdown; /** * The up quark mass */ Energy _mup; /** * The strange quark mass */ Energy _mstrange; /** * The charm quark mass */ Energy _mcharm; /** * The bottom quark mass */ Energy _mbottom; /** * The masses of the quarks as a vector */ vector _mquark; //@} /** @name Wave function parameters */ //@{ /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{d}\f$ */ Energy _betaSud; /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{s}\f$ */ Energy _betaSus; /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{c}\f$ */ Energy _betaSuc; /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{b}\f$ */ Energy _betaSub; /** * The s-wave variational parameters as a vector. */ vector > _betaS; /** * The wavefunction p-wave \f$\beta\f$ variational parameters for \f$u\bar{d}\f$ */ Energy _betaPud; /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{s}\f$ */ Energy _betaPus; /** * The wavefunction s-wave \f$\beta\f$ variational parameters for \f$u\bar{c}\f$ */ Energy _betaPuc; /** * The p-wave variational parameters as a vector */ vector > _betaP; //@} /** * The \f$\eta-\eta'\f$ mixing angle */ double _thetaeta; }; } #endif /* HERWIG_ISGWFormFactor_H */ 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/SMWWHVertex.h b/Models/StandardModel/SMWWHVertex.h --- a/Models/StandardModel/SMWWHVertex.h +++ b/Models/StandardModel/SMWWHVertex.h @@ -1,135 +1,136 @@ // -*- C++ -*- // // SMWWHVertex.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_SMWWHVertex_H #define HERWIG_SMWWHVertex_H // // This is the declaration of the SMWWHVertex class. // #include "ThePEG/Helicity/Vertex/Scalar/VVSVertex.h" #include "ThePEG/PDT/EnumParticles.h" #include "Herwig/Models/StandardModel/StandardModel.h" namespace Herwig { using namespace ThePEG; /** \ingroup Helicity * * The SMWWHVertex is the implementation of the * coupling of two electroweak gauge bosons to the Higgs in the Standard * Model. It inherits from VVSVertex and implements the setCoupling member. * * @see VVSVertex * @see VertexBase */ class SMWWHVertex: public VVSVertex { public: /** * Default constructor. */ SMWWHVertex(); /** * Calculate the couplings. * @param q2 The scale \f$q^2\f$ for the coupling at the vertex. * @param part1 The ParticleData pointer for the first particle. * @param part2 The ParticleData pointer for the second particle. * @param part3 The ParticleData pointer for the third particle. */ virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3); public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * Standard Init function used to initialize the interfaces. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const {return new_ptr(*this);} /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const {return new_ptr(*this);} //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving and * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); //@} private: /** + * Private and non-existent assignment operator. */ SMWWHVertex & operator=(const SMWWHVertex &) = delete; private: /** * Storage of the couplings. */ //@{ /** * The last value of the electroweak coupling calculated. */ Complex _couplast; /** * The scale \f$q^2\f$ at which the coupling was last evaluated. */ Energy2 _q2last; /** * The mass of the \f$W\f$ boson. */ Energy _mw; /** * The factor for the \f$Z\f$ vertex. */ double _zfact; //@} }; } #endif /* HERWIG_SMWWHVertex_H */ diff --git a/Models/Transplanckian/METRP2to2.h b/Models/Transplanckian/METRP2to2.h --- a/Models/Transplanckian/METRP2to2.h +++ b/Models/Transplanckian/METRP2to2.h @@ -1,235 +1,237 @@ // -*- C++ -*- // // METRP2to2.h is a part of Herwig - A multi-purpose Monte Carlo event generator // Copyright (C) 2009-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_METRP2to2_H #define HERWIG_METRP2to2_H // // This is the declaration of the METRP2to2 class. // // The implementation of this process is based upon hep-ph/0112161 by G.F. Giudice, R. Rattazzi, J.D. Wells. #include "Herwig/MatrixElement/HwMEBase.h" #include "ThePEG/Repository/UseRandom.h" #include "Herwig/Utilities/Interpolator.h" namespace Herwig { using namespace ThePEG; /** * The METRP2to2 class implements the matrix elements for * Transplanckian \f$2\to2\f$ scattering process + * + * @see \ref METRP2to2Interfaces "The interfaces" defined for METRP2to2. */ class METRP2to2: public HwMEBase { public: /** * The default constructor. */ METRP2to2(); /** @name Virtual functions required by the MEBase class. */ //@{ /** * Return the order in \f$\alpha_S\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaS() const { return 0; } /** * Return the order in \f$\alpha_{EW}\f$ in which this matrix * element is given. */ virtual unsigned int orderInAlphaEW() const { return 0; } /** * Return the scale associated with the last set phase space point. */ virtual Energy2 scale() const; /** * Add all possible diagrams with the add() function. */ virtual void getDiagrams() const; /** * Get diagram selector. With the information previously supplied with the * setKinematics method, a derived class may optionally * override this method to weight the given diagrams with their * (although certainly not physical) relative probabilities. * @param dv the diagrams to be weighted. * @return a Selector relating the given diagrams to their weights. */ virtual Selector diagrams(const DiagramVector & dv) const; /** * Return a Selector with possible colour geometries for the selected * diagram weighted by their relative probabilities. * @param diag the diagram chosen. * @return the possible colour geometries weighted by their * relative probabilities. */ virtual Selector colourGeometries(tcDiagPtr diag) const; /** * The matrix element for the kinematical configuration * previously provided by the last call to setKinematics(), suitably * scaled by sHat() to give a dimension-less number. * @return the matrix element scaled with sHat() to give a * dimensionless number. */ virtual double me2() const; //@} protected: /** @name Standard Interfaced functions. */ //@{ /** * Initialize this object after the setup phase before saving an * EventGenerator to disk. * @throws InitException if object could not be initialized properly. */ virtual void doinit(); /** * Rebind pointer to other Interfaced objects. Called in the setup phase * after all objects used in an EventGenerator has been cloned so that * the pointers will refer to the cloned objects afterwards. * @param trans a TranslationMap relating the original objects to * their respective clones. * @throws RebindException if no cloned object was found for a given * pointer. */ virtual void rebind(const TranslationMap & trans); /** * Return a vector of all pointers to Interfaced objects used in this * object. * @return a vector of pointers. */ virtual IVector getReferences(); //@} public: /** @name Functions used by the persistent I/O system. */ //@{ /** * Function used to write out object persistently. * @param os the persistent output stream written to. */ void persistentOutput(PersistentOStream & os) const; /** * Function used to read in object persistently. * @param is the persistent input stream read from. * @param version the version number of the object when written. */ void persistentInput(PersistentIStream & is, int version); //@} /** * The standard Init function used to initialize the interfaces. * Called exactly once for each class by the class description system * before the main function starts or * when this class is dynamically loaded. */ static void Init(); protected: /** @name Clone Methods. */ //@{ /** * Make a simple clone of this object. * @return a pointer to the new object. */ virtual IBPtr clone() const; /** Make a clone of this object, possibly modifying the cloned object * to make it sane. * @return a pointer to the new object. */ virtual IBPtr fullclone() const; //@} protected: /** @name Helper functions for me2. */ //@{ /** * The function which calculates b_c according to hep-ph/0112161, eq.(18) */ InvEnergy bccalc(Energy2 s) const; /** * A_ny calculates part of the matrix element squared (divided by 16 * pi^2) according to hep-ph/0112161. eq.(19) */ double A_ny(Energy2 s, Energy2 t) const; /** * fpoint initializes the matrix of pre-calculated points for the functions F_n(y) (hep-ph/0112161, eq.(20) */ double fpoint(double x) const; /** * The asymptotic form of the F_n(y) functions, used for y>20, according to hep-ph/0112161, eq. (25) */ double fnyasympt(double y) const; //@} /** * Setup the interpolator */ void setup_interpolator(); private: /** * The assignment operator is private and must never be called. * In fact, it should not even be implemented. */ METRP2to2 & operator=(const METRP2to2 &) = delete; private: /** * Interpolator */ Interpolator::Ptr _interpol; /** * Maximum number of quark flavours to include */ unsigned int _maxflavour; /** * Number of Extra dimensions (>=2) */ unsigned int _ndim; /** * The Extra-dimensional Planck mass */ Energy _planckmass; /** * Processes to include */ unsigned int _process; }; } #endif /* HERWIG_METRP2to2_H */ diff --git a/Shower/Dipole/Kernels/FIgx2ggxDipoleKernel.cc b/Shower/Dipole/Kernels/FIgx2ggxDipoleKernel.cc --- a/Shower/Dipole/Kernels/FIgx2ggxDipoleKernel.cc +++ b/Shower/Dipole/Kernels/FIgx2ggxDipoleKernel.cc @@ -1,116 +1,116 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the FIgx2ggxDipoleKernel class. // #include "FIgx2ggxDipoleKernel.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" using namespace Herwig; FIgx2ggxDipoleKernel::FIgx2ggxDipoleKernel() : DipoleSplittingKernel() {} FIgx2ggxDipoleKernel::~FIgx2ggxDipoleKernel() {} IBPtr FIgx2ggxDipoleKernel::clone() const { return new_ptr(*this); } IBPtr FIgx2ggxDipoleKernel::fullclone() const { return new_ptr(*this); } bool FIgx2ggxDipoleKernel::canHandle(const DipoleIndex& ind) const { return useThisKernel() && ind.emitterData()->id() == ParticleID::g && ind.spectatorData()->mass() == ZERO && !ind.initialStateEmitter() && ind.initialStateSpectator(); } bool FIgx2ggxDipoleKernel::canHandleEquivalent(const DipoleIndex& a, const DipoleSplittingKernel& sk, const DipoleIndex& b) const { assert(canHandle(a)); if ( !canHandle(b) ) return false; return sk.emitter(b)->id() == ParticleID::g && sk.emission(b)->id() == ParticleID::g && a.spectatorPDF() == b.spectatorPDF(); } tcPDPtr FIgx2ggxDipoleKernel::emitter(const DipoleIndex&) const { return getParticleData(ParticleID::g); } tcPDPtr FIgx2ggxDipoleKernel::emission(const DipoleIndex&) const { return getParticleData(ParticleID::g); } tcPDPtr FIgx2ggxDipoleKernel::spectator(const DipoleIndex& ind) const { return ind.spectatorData(); } double FIgx2ggxDipoleKernel::evaluate(const DipoleSplittingInfo& split) const { double ret = alphaPDF(split); double z = split.lastZ(); double x = 1. / ( 1. + sqr(split.lastPt()/split.scale()) / (z*(1.-z)) ); - double S1=1./(1.-z*(1.-x)); - double S2=1./(1.-(1.-z)*(1.-x)); + double S1=1./(1.-z+(1.-x)); + double S2=1./(1.-(1.-z)+(1.-x)); double NS=(-2 + z*(1.-z)+(1.-x)*(1.+x*z*(1.-z))); if( theAsymmetryOption == 0 ){ ret *= 3.*( S1 + 0.5 * NS); }else if ( theAsymmetryOption == 1 ){ ret *= 3.*z*( S1 +S2 + NS ); }else{ ret *= 3.*0.5*( S1 + S2 + NS ); } return ret > 0. ? ret : 0.; } // If needed, insert default implementations of function defined // in the InterfacedBase class here (using ThePEG-interfaced-impl in Emacs). void FIgx2ggxDipoleKernel::persistentOutput(PersistentOStream & os) const { os << theAsymmetryOption; } void FIgx2ggxDipoleKernel::persistentInput(PersistentIStream & is, int) { is >> theAsymmetryOption; } ClassDescription FIgx2ggxDipoleKernel::initFIgx2ggxDipoleKernel; // Definition of the static class description member. void FIgx2ggxDipoleKernel::Init() { static ClassDocumentation documentation ("FIgx2ggxDipoleKernel"); static Parameter interfacetheAsymmetryOption ("AsymmetryOption", "The asymmetry option for final state gluon spliitings.", &FIgx2ggxDipoleKernel::theAsymmetryOption, 1, 0, 0, false, false, Interface::lowerlim); } diff --git a/src/defaults/HerwigCleanup.in b/src/defaults/HerwigCleanup.in deleted file mode 100644 --- a/src/defaults/HerwigCleanup.in +++ /dev/null @@ -1,4 +0,0 @@ -# -*- ThePEG-repository -*- - -# Clean up all Herwig related directories -rrmdir /Herwig diff --git a/src/defaults/Makefile.am b/src/defaults/Makefile.am --- a/src/defaults/Makefile.am +++ b/src/defaults/Makefile.am @@ -1,47 +1,47 @@ BUILT_SOURCES = done-all-links defaultsdir = ${pkgdatadir}/defaults INPUTFILES = Analysis.in \ baryon_decays.in baryons.in boson_decays.in \ bosons.in Cuts.in decayers.in \ Decays.in DiffractiveParticles.in diquarks.in \ -Hadronization.in HerwigDefaults.in HerwigCleanup.in \ +Hadronization.in HerwigDefaults.in \ lepton_decays.in leptons.in \ masses.in MatrixElements.in meson_decays.in mesons.in Model.in BSM.in \ Particles.in QEDRadiation.in quark_decays.in quarks.in \ setup.gosam.in Shower.in StandardModelVertices.in UnderlyingEvent.in widths.in Partons.in \ MatchboxDefaults.in DipoleShowerDefaults.in MatchboxLoopInduced.in Samplers.in \ EvtGenDecayer.in EvtGenBDecays.in HerwigBDecays.in MatchboxMergingDefaults.in dist_defaults_DATA = $(INPUTFILES) EXTRA_DIST = PDF.in.in defaults_DATA = PDF.in CLEANFILES = PDF.in done-all-links ## For an explanation of this magic, see autoconf book 4.7.2 edit = sed -e "s,@HERWIG_PDF_POMERON\@,`cd $(top_srcdir) && pwd`/PDF/diffraction/," ## may still be broken and need $(DESTDIR) installnamePOMERON = $(pkgdatadir)/PDF/diffraction/ install-data-hook: rm -f $(DESTDIR)$(defaultsdir)/PDF.in sed -e 's,@HERWIG_PDF_POMERON\@,$(installnamePOMERON),' \ $(srcdir)/PDF.in.in > $(DESTDIR)$(defaultsdir)/PDF.in PDF.in: Makefile $(srcdir)/PDF.in.in @echo "Updating PDF.in" @rm -f PDF.in PDF.in.tmp @$(edit) $(srcdir)/PDF.in.in > PDF.in.tmp @mv PDF.in.tmp PDF.in done-all-links: $(INPUTFILES) @echo "Linking input files" @for i in $(INPUTFILES); do \ if test -f $(srcdir)/$$i -a ! -e $$i; then \ $(LN_S) -f $(srcdir)/$$i; fi; done @touch done-all-links