diff --git a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc --- a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc +++ b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc @@ -1,360 +1,360 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the VectorMeson2BaryonsDecayer class. // #include "VectorMeson2BaryonsDecayer.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/ParVector.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/PDT/DecayMode.h" #include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/Decay/TwoBodyDecayMatrixElement.h" #include "ThePEG/Helicity/HelicityFunctions.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/RSSpinorBarWaveFunction.h" using namespace Herwig; using namespace ThePEG::Helicity; void VectorMeson2BaryonsDecayer::doinitrun() { DecayIntegrator::doinitrun(); if(initialize()) { for(unsigned int ix=0;ixmaxWeight(); } } } void VectorMeson2BaryonsDecayer::doinit() { DecayIntegrator::doinit(); // check the parameters arew consistent unsigned int isize=gm_.size(); if(isize!=incoming_.size() || isize!=outgoingf_.size()|| isize!=outgoinga_.size() || isize!= maxweight_.size()|| isize!= phi_.size() || isize!= ge_.size()) throw InitException() << "Inconsistent parameters in VectorMeson2" << "BaryonsDecayer::doiin() " << Exception::runerror; // set up the integration channels PhaseSpaceModePtr mode; for(unsigned int ix=0;ixid()); int idbar = parent->CC() ? parent->CC()->id() : id; int id1(children[0]->id()); int id1bar = children[0]->CC() ? children[0]->CC()->id() : id1; int id2(children[1]->id()); int id2bar = children[1]->CC() ? children[1]->CC()->id() : id2; int imode(-1); unsigned int ix(0); cc=false; do { if(incoming_[ix]==id ) { if((id1 ==outgoingf_[ix]&&id2 ==outgoinga_[ix])|| (id2 ==outgoingf_[ix]&&id1 ==outgoinga_[ix])) imode=ix; } if(incoming_[ix]==idbar) { if((id1bar==outgoingf_[ix]&&id2bar==outgoinga_[ix])|| (id2bar==outgoingf_[ix]&&id1bar==outgoinga_[ix])) { imode=ix; cc=true; } } ++ix; } while(imode<0&&ix> gm_ >> ge_ >> phi_ >> incoming_ >> outgoingf_ >> outgoinga_ >> maxweight_; } DescribeClass describeHerwigVectorMeson2BaryonsDecayer("Herwig::VectorMeson2BaryonsDecayer", "HwVMDecay.so"); void VectorMeson2BaryonsDecayer::Init() { static ClassDocumentation documentation ("The VectorMeson2BaryonsDecayer class is designed for " "the decay of vector mesons to baryons."); static ParVector interfaceIncoming ("Incoming", "The PDG code for the incoming particle", &VectorMeson2BaryonsDecayer::incoming_, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcoming1 ("OutgoingBaryons", "The PDG code for the outgoing fermion", &VectorMeson2BaryonsDecayer::outgoingf_, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceOutcoming2 ("OutgoingAntiBaryons", "The PDG code for the second outgoing anti-fermion", &VectorMeson2BaryonsDecayer::outgoinga_, 0, 0, 0, -10000000, 10000000, false, false, true); static ParVector interfaceGM ("GM", "The value of the GM form factor", &VectorMeson2BaryonsDecayer::gm_, -1, 0., -1000., 1000., false, false, Interface::limited); static ParVector interfaceGE ("GE", "The value of the GE form factor", &VectorMeson2BaryonsDecayer::ge_, -1, 0., -1000., 1000., false, false, Interface::limited); static ParVector interfacePhi ("Phi", "The phase of the GE form factor", &VectorMeson2BaryonsDecayer::phi_, -1, 0., -Constants::pi, Constants::pi, false, false, Interface::limited); static ParVector interfaceMaxWeight ("MaxWeight", "The maximum weight for the decay mode", &VectorMeson2BaryonsDecayer::maxweight_, 0, 0, 0, -10000000, 10000000, false, false, true); } void VectorMeson2BaryonsDecayer:: constructSpinInfo(const Particle & part, ParticleVector decay) const { unsigned int iferm(0),ianti(1); if(outgoingf_[imode()]!=decay[iferm]->id()) swap(iferm,ianti); VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast(&part), incoming,true,false); // outgoing fermion if(decay[iferm]->dataPtr()->iSpin()==PDT::Spin1Half) SpinorBarWaveFunction:: constructSpinInfo(wavebar_,decay[iferm],outgoing,true); else RSSpinorBarWaveFunction:: constructSpinInfo(wave2bar_,decay[iferm],outgoing,true); // outgoing antifermion if(decay[ianti]->dataPtr()->iSpin()==PDT::Spin1Half) SpinorWaveFunction:: constructSpinInfo(wave_ ,decay[ianti],outgoing,true); else RSSpinorWaveFunction:: constructSpinInfo(wave2_,decay[ianti],outgoing,true); } double VectorMeson2BaryonsDecayer::me2(const int,const Particle & part, const tPDVector & outgoing, const vector & momenta, MEOption meopt) const { // initialze me if(!ME()) ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,outgoing[0]->iSpin(),outgoing[0]->iSpin()))); // fermion and antifermion unsigned int iferm(0),ianti(1); if(outgoingf_[imode()]!=outgoing[iferm]->id()) swap(iferm,ianti); // initialization if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_,rho_, const_ptr_cast(&part), incoming,false); } // spin 1/2 if(outgoing[0]->iSpin()==PDT::Spin1Half && outgoing[1]->iSpin()==PDT::Spin1Half) { wave_.resize(2); wavebar_.resize(2); for(unsigned int ix=0;ix<2;++ix) { wavebar_[ix] = HelicityFunctions::dimensionedSpinorBar(-momenta[iferm],ix,Helicity::outgoing); wave_ [ix] = HelicityFunctions::dimensionedSpinor (-momenta[ianti],ix,Helicity::outgoing); } // coefficients Complex GM = gm_[imode()]; Complex GE = ge_[imode()]*exp(Complex(0.,phi_[imode()])); LorentzPolarizationVector c2 = -2.*outgoing[0]->mass()/(4.*sqr(outgoing[0]->mass())-sqr(part.mass()))* (GM-GE)*(momenta[iferm]-momenta[ianti]); // now compute the currents LorentzPolarizationVector temp; //double mesum(0.); for(unsigned ix=0;ix<2;++ix) { for(unsigned iy=0;iy<2;++iy) { LorentzPolarizationVector temp = (GM*wave_[ix].vectorCurrent(wavebar_[iy])+c2*wave_[ix].scalar(wavebar_[iy]))/part.mass(); for(unsigned int iz=0;iz<3;++iz) { if(iferm>ianti) (*ME())(iz,ix,iy)=vectors_[iz].dot(temp); else (*ME())(iz,iy,ix)=vectors_[iz].dot(temp); //mesum += norm(vectors_[iz].dot(temp)); } } } double me = ME()->contract(rho_).real(); // cerr << "testing decay " << part.PDGName() << " -> " << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << "\n"; // cerr << "testing ME " << mesum/3. << " " << me << " " << 4./3.*(norm(GM)+2.*sqr(outgoing[0]->mass()/part.mass())*norm(GE)) << "\n"; // cerr << "testing gamma " << mesum/3./8./Constants::pi*sqrt(sqr(part.mass())-4.*sqr(outgoing[0]->mass()))/MeV << "\n"; // return the answer return me; } // spin 3/2 else if(outgoing[0]->iSpin()==PDT::Spin3Half && outgoing[0]->iSpin()==PDT::Spin3Half) { wave2_.resize(4); wave2bar_.resize(4); wave_.resize(4); wavebar_.resize(4); RSSpinorBarWaveFunction swave(momenta[iferm],outgoing[iferm],Helicity::outgoing); RSSpinorWaveFunction awave(momenta[ianti],outgoing[ianti],Helicity::outgoing); LorentzPolarizationVector vtemp = part.momentum()/part.mass(); for(unsigned int ix=0;ix<4;++ix) { swave.reset(ix); awave.reset(ix); wave2bar_[ix] = swave.dimensionedWf(); wavebar_ [ix] = wave2bar_[ix].dot(vtemp); wave2_ [ix] = awave.dimensionedWf(); wave_ [ix] = wave2_[ix].dot(vtemp); } // coefficients Complex GM = gm_[imode()]; Complex GE = ge_[imode()]*exp(Complex(0.,phi_[imode()])); LorentzPolarizationVector c2 = -2.*outgoing[0]->mass()/(4.*sqr(outgoing[0]->mass())-sqr(part.mass()))* (GM-GE)*(momenta[iferm]-momenta[ianti]); // now compute the currents for(unsigned ix=0;ix<4;++ix) { for(unsigned iy=0;iy<4;++iy) { // q(al)q(be) piece LorentzPolarizationVector temp2 = (GM*wave_[ix].vectorCurrent(wavebar_[iy])+c2*wave_[ix].scalar(wavebar_[iy]))* 2.*part.mass()/(4.*sqr(outgoing[0]->mass())-sqr(part.mass())); // g(al)g(be) * GM-GE piece LorentzPolarizationVector temp3 = wave2_[ix].generalScalar(wave2bar_[iy],1.,1.)*c2/part.mass(); // g(al)g(be) * gamma^mu LorentzPolarizationVector temp1(GM/part.mass()*(wave2bar_[iy](0,3)*wave2_[ix](0,0) + wave2bar_[iy](0,2)*wave2_[ix](0,1) - wave2bar_[iy](0,1)*wave2_[ix](0,2) - wave2bar_[iy](0,0)*wave2_[ix](0,3) + wave2bar_[iy](1,3)*wave2_[ix](1,0) + wave2bar_[iy](1,2)*wave2_[ix](1,1) - wave2bar_[iy](1,1)*wave2_[ix](1,2) - wave2bar_[iy](1,0)*wave2_[ix](1,3) + wave2bar_[iy](2,3)*wave2_[ix](2,0) + wave2bar_[iy](2,2)*wave2_[ix](2,1) - wave2bar_[iy](2,1)*wave2_[ix](2,2) - wave2bar_[iy](2,0)*wave2_[ix](2,3) - wave2bar_[iy](3,3)*wave2_[ix](3,0) - wave2bar_[iy](3,2)*wave2_[ix](3,1) + wave2bar_[iy](3,1)*wave2_[ix](3,2) + wave2bar_[iy](3,0)*wave2_[ix](3,3)), Complex(0,1)*GM/part.mass()*(wave2bar_[iy](0,3)*wave2_[ix](0,0) - wave2bar_[iy](0,2)*wave2_[ix](0,1) - wave2bar_[iy](0,1)*wave2_[ix](0,2) + wave2bar_[iy](0,0)*wave2_[ix](0,3) + wave2bar_[iy](1,3)*wave2_[ix](1,0) - wave2bar_[iy](1,2)*wave2_[ix](1,1) - wave2bar_[iy](1,1)*wave2_[ix](1,2) + wave2bar_[iy](1,0)*wave2_[ix](1,3) + wave2bar_[iy](2,3)*wave2_[ix](2,0) - wave2bar_[iy](2,2)*wave2_[ix](2,1) - wave2bar_[iy](2,1)*wave2_[ix](2,2) + wave2bar_[iy](2,0)*wave2_[ix](2,3) - wave2bar_[iy](3,3)*wave2_[ix](3,0) + wave2bar_[iy](3,2)*wave2_[ix](3,1) + wave2bar_[iy](3,1)*wave2_[ix](3,2) - wave2bar_[iy](3,0)*wave2_[ix](3,3)), GM/part.mass()*(wave2bar_[iy](0,2)*wave2_[ix](0,0) - wave2bar_[iy](0,3)*wave2_[ix](0,1) - wave2bar_[iy](0,0)*wave2_[ix](0,2) + wave2bar_[iy](0,1)*wave2_[ix](0,3) + wave2bar_[iy](1,2)*wave2_[ix](1,0) - wave2bar_[iy](1,3)*wave2_[ix](1,1) - wave2bar_[iy](1,0)*wave2_[ix](1,2) + wave2bar_[iy](1,1)*wave2_[ix](1,3) + wave2bar_[iy](2,2)*wave2_[ix](2,0) - wave2bar_[iy](2,3)*wave2_[ix](2,1) - wave2bar_[iy](2,0)*wave2_[ix](2,2) + wave2bar_[iy](2,1)*wave2_[ix](2,3) - wave2bar_[iy](3,2)*wave2_[ix](3,0) + wave2bar_[iy](3,3)*wave2_[ix](3,1) + wave2bar_[iy](3,0)*wave2_[ix](3,2) - wave2bar_[iy](3,1)*wave2_[ix](3,3)), GM/part.mass()*(-wave2bar_[iy](0,2)*wave2_[ix](0,0) - wave2bar_[iy](0,3)*wave2_[ix](0,1) - wave2bar_[iy](0,0)*wave2_[ix](0,2) - wave2bar_[iy](0,1)*wave2_[ix](0,3) - wave2bar_[iy](1,2)*wave2_[ix](1,0) - wave2bar_[iy](1,3)*wave2_[ix](1,1) - wave2bar_[iy](1,0)*wave2_[ix](1,2) - wave2bar_[iy](1,1)*wave2_[ix](1,3) - wave2bar_[iy](2,2)*wave2_[ix](2,0) - wave2bar_[iy](2,3)*wave2_[ix](2,1) - wave2bar_[iy](2,0)*wave2_[ix](2,2) - wave2bar_[iy](2,1)*wave2_[ix](2,3) + wave2bar_[iy](3,2)*wave2_[ix](3,0) + wave2bar_[iy](3,3)*wave2_[ix](3,1) + wave2bar_[iy](3,0)*wave2_[ix](3,2) + wave2bar_[iy](3,1)*wave2_[ix](3,3))); LorentzPolarizationVector temp = temp1+temp2+temp3; for(unsigned int iz=0;iz<3;++iz) { if(iferm>ianti) (*ME())(iz,ix,iy)=vectors_[iz].dot(temp); else (*ME())(iz,iy,ix)=vectors_[iz].dot(temp); } } } // double mesum = ME()->contract(RhoDMatrix(PDT::Spin1)).real(); // generator()->log() << "testing decay " << part.PDGName() << " -> " << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << "\n"; // generator()->log() << "testing ME " << mesum << " " << me << " " << 1./3.*(40.*norm(GM)/9.+16.*sqr(outgoing[0]->mass()/part.mass())*norm(GE)) << "\n"; // return the answer return ME()->contract(rho_).real(); } else assert(false); } // output the setup information for the particle database void VectorMeson2BaryonsDecayer::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DecayIntegrator base class DecayIntegrator::dataBaseOutput(output,false); // the rest of the parameters for(unsigned int ix=0;ix