diff --git a/Decay/Dalitz/ScalarTo3ScalarDalitz.cc b/Decay/Dalitz/ScalarTo3ScalarDalitz.cc --- a/Decay/Dalitz/ScalarTo3ScalarDalitz.cc +++ b/Decay/Dalitz/ScalarTo3ScalarDalitz.cc @@ -1,178 +1,174 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the ScalarTo3ScalarDalitz class. // #include "ScalarTo3ScalarDalitz.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/Interface/Switch.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; using namespace ThePEG::Helicity; IBPtr ScalarTo3ScalarDalitz::clone() const { return new_ptr(*this); } IBPtr ScalarTo3ScalarDalitz::fullclone() const { return new_ptr(*this); } void ScalarTo3ScalarDalitz::persistentOutput(PersistentOStream & os) const { os << useResonanceMass_ ; } void ScalarTo3ScalarDalitz::persistentInput(PersistentIStream & is, int) { is >> useResonanceMass_; } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigScalarTo3ScalarDalitz("Herwig::ScalarTo3ScalarDalitz", "HwDalitzDecay.so"); void ScalarTo3ScalarDalitz::Init() { static ClassDocumentation documentation ("The ScalarTo3ScalarDalitz class provides a base class for " "weak three-body decays of bottom and charm mesons"); static Switch interfaceResonanceMass ("ResonanceMass", "Whether to use the kinematic mass or the resonance pole mass for the denominator in kinematic expressions", &ScalarTo3ScalarDalitz::useResonanceMass_, false, false, false); static SwitchOption interfaceResonanceMassYes (interfaceResonanceMass, "Yes", "Use the resonance mass, to be avoided only use if do in experimental fit", true); static SwitchOption interfaceResonanceMassNo (interfaceResonanceMass, "No", "Use the correct kinematic mass", false); } void ScalarTo3ScalarDalitz:: constructSpinInfo(const Particle & part, ParticleVector decay) const { // set up the spin information for the decay products ScalarWaveFunction::constructSpinInfo(const_ptr_cast(&part), incoming,true); for(unsigned int ix=0;ix<3;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); } double ScalarTo3ScalarDalitz::me2(const int ichan, const Particle & part, const tPDVector & , const vector & momenta, MEOption meopt) const { static const Complex ii(0.,1.); if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin0,PDT::Spin0))); useMe(); if(meopt==Initialize) { ScalarWaveFunction:: calculateWaveFunctions(rho_,const_ptr_cast(&part),incoming); } // set the kinematics mD_ = part.mass(); for(unsigned int ix=0;ixamp; if (resonances()[i]->type==ResonanceType::NonResonant) return output; // mass of the resonance const Energy & mR = resonances()[i]->mass ; // locations of the outgoing particles const unsigned int &d1 = resonances()[i]->daughter1; const unsigned int &d2 = resonances()[i]->daughter2; const unsigned int &sp = resonances()[i]->spectator; // compute the Breit-Wigner times resonance form factor piece output *= resonances()[i]->BreitWigner(m2_[d1][d2],mOut_[d1],mOut_[d2]); // angular piece // mass for the denominator Energy mDen = useResonanceMass_ ? resonances()[i]->mass : m2_[d1][d2]; // denominator for the older form of the amplitude Energy2 denom = GeV2; if (resonances()[i]->type/10 == 1 ) { Energy2 pa2 = 0.25*(sqr(m2_[d1][d2])-2.*(sqr(mOut_[d1])+sqr(mOut_[d2])) + sqr(sqr(mOut_[d1])-sqr(mOut_[d2]))/sqr(m2_[d1][d2])); Energy2 pc2 = 0.25*(sqr(m2_[d1][d2])-2.*(sqr(mD_ )+sqr(mOut_[sp])) + sqr(sqr(mD_ )-sqr(mOut_[sp]))/sqr(m2_[d1][d2])); denom = 4.*sqrt(pa2*pc2); } // vectors if(abs(resonances()[i]->type)%10==3) { output *= (sqr(m2_[d2][sp])-sqr(m2_[d1][sp]) + ( sqr(mD_)-sqr(mOut_[sp]))*(sqr(mOut_[d1])-sqr(mOut_[d2]))/sqr(mDen) )/denom; } else if(abs(resonances()[i]->type)%10==5) { output *= 1./sqr(denom)*( sqr( sqr(m2_[d2][sp]) - sqr(m2_[d1][sp]) +(sqr(mD_)-sqr(mOut_[sp]))*(sqr(mOut_[d1])-sqr(mOut_[d2]))/(sqr(mDen))) - (sqr(m2_[d1][d2])-2* sqr(mD_)-2*sqr(mOut_[sp]) + sqr((sqr( mD_) - sqr(mOut_[sp]))/mDen))* (sqr(m2_[d1][d2])-2*sqr(mOut_[d1])-2*sqr(mOut_[d2]) + sqr((sqr(mOut_[d1]) - sqr(mOut_[d2]))/mDen))/3.); } // spin zero and non-resonant is done now if((abs(resonances()[i]->type)%10==1 && resonances()[i]->type != ResonanceType::Spin0Gauss ) || abs(resonances()[i]->type/10)==10) { return output; } // spin piece x Blatt-Weisskopf for parent else { double fD=1.; // for the D decay Energy pD = sqrt(max(ZERO,(0.25*sqr(sqr(mD_)-sqr(mR)-sqr(mOut_[sp])) - sqr(mR*mOut_[sp]))/sqr(mD_))); Energy pDAB= sqrt( 0.25*sqr(sqr(mD_)-sqr(m2_[d1][d2])-sqr(mOut_[sp])) - sqr(m2_[d1][d2]*mOut_[sp]))/mD_; double r2A(parentRadius() *pD),r2B(parentRadius() *pDAB); // Blatt-Weisskopf factors and spin piece switch (resonances()[i]->type) { case ResonanceType::Spin0Gauss: fD = exp(-(r2B-r2A)/12.); - // cerr << "testing scalar " << i << " " << mR/GeV << " " << exp(r2A/12.) << "\n"; - output *= fD; break; case ResonanceType::Spin1: case ResonanceType::Spin1E691 : case ResonanceType::Spin1GS : fD=sqrt( (1. + sqr(r2A)) / (1. + sqr(r2B)) ); - // cerr << "testing vector " << i << " " << mR/GeV << " " << sqrt(1. + sqr(r2A)) << "\n"; break; case ResonanceType::Spin2: case ResonanceType::Spin2E691: fD = sqrt( (9. + sqr(r2A)*(3.+sqr(r2A))) / (9. + sqr(r2B)*(3.+sqr(r2B)))); - // cerr << "testing tensor " << i << " " << mR/GeV << " " << sqrt( (9. + sqr(r2A)*(3.+sqr(r2A)))) << "\n"; break; default : assert(false); } output *= fD; return output; } } void ScalarTo3ScalarDalitz::dataBaseOutput(ofstream & output, bool header) const { if(header) output << "update decayers set parameters=\""; // parameters for the DalitzBase base class DalitzBase::dataBaseOutput(output,false); output << "newdef " << name() << ":ResonanceMass " << useResonanceMass_ << "\n"; output << "\n"; if(header) { output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl; } } diff --git a/Decay/Dalitz/VectorTo3PseudoScalarDalitz.cc b/Decay/Dalitz/VectorTo3PseudoScalarDalitz.cc --- a/Decay/Dalitz/VectorTo3PseudoScalarDalitz.cc +++ b/Decay/Dalitz/VectorTo3PseudoScalarDalitz.cc @@ -1,221 +1,206 @@ // -*- C++ -*- // // This is the implementation of the non-inlined, non-templated member // functions of the VectorTo3PseudoScalarDalitz class. // #include "VectorTo3PseudoScalarDalitz.h" #include "ThePEG/Interface/ClassDocumentation.h" #include "ThePEG/Interface/Parameter.h" #include "ThePEG/EventRecord/Particle.h" #include "ThePEG/Repository/UseRandom.h" #include "ThePEG/Repository/EventGenerator.h" #include "ThePEG/Utilities/DescribeClass.h" #include "ThePEG/Persistency/PersistentOStream.h" #include "ThePEG/Persistency/PersistentIStream.h" #include "ThePEG/Helicity/epsilon.h" #include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h" #include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h" #include "Herwig/PDT/ThreeBodyAllOnCalculator.h" #include "Herwig/Decay/GeneralDecayMatrixElement.h" using namespace Herwig; IBPtr VectorTo3PseudoScalarDalitz::clone() const { return new_ptr(*this); } IBPtr VectorTo3PseudoScalarDalitz::fullclone() const { return new_ptr(*this); } void VectorTo3PseudoScalarDalitz::persistentOutput(PersistentOStream & os) const { os << ounit(coupling_,1./GeV); } void VectorTo3PseudoScalarDalitz::persistentInput(PersistentIStream & is, int) { is >> iunit(coupling_,1./GeV); } // The following static variable is needed for the type // description system in ThePEG. DescribeClass describeHerwigVectorTo3PseudoScalarDalitz("Herwig::VectorTo3PseudoScalarDalitz", "HwDalitzDecay.so"); void VectorTo3PseudoScalarDalitz::Init() { static ClassDocumentation documentation ("The VectorTo3PseudoScalarDalitz class provides a base class " "for the decay of vector mesons to 3 pseudoscalar mesons"); static Parameter interfaceCouopling ("Coupling", "The coupling for the normalisation of the mode", &VectorTo3PseudoScalarDalitz::coupling_, 1./GeV, 1./GeV, 0./GeV, 1000./GeV, false, false, Interface::limited); } void VectorTo3PseudoScalarDalitz:: constructSpinInfo(const Particle & part, ParticleVector decay) const { VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast(&part), incoming,true,false); // set up the spin information for the decay products for(unsigned int ix=0;ix<3;++ix) ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true); } double VectorTo3PseudoScalarDalitz::me2(const int ichan, const Particle & part, const tPDVector & , const vector & momenta, MEOption meopt) const { if(!ME()) ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0,PDT::Spin0))); useMe(); if(meopt==Initialize) { VectorWaveFunction::calculateWaveFunctions(vectors_,rho_, const_ptr_cast(&part), incoming,false); } // set the kinematics mD_ = part.mass(); for(unsigned int ix=0;ix amp = amplitude(ichan); // polarization vector piece LorentzVector > scalar = epsilon(momenta[0],momenta[1],momenta[2]); // compute the matrix element for(unsigned int ix=0;ix<3;++ix) { (*ME())(ix,0,0,0) = Complex(coupling_*amp*scalar.dot(vectors_[ix])); } // return the answer return (ME()->contract(rho_)).real(); } double VectorTo3PseudoScalarDalitz:: threeBodyMatrixElement(const int , const Energy2 q2, const Energy2 s3, const Energy2 s2, const Energy2 s1, const Energy m1, const Energy m2, const Energy m3) const { mD_ = sqrt(q2); mOut_[0] = m1; mOut_[1] = m2; mOut_[2] = m3; m2_[0][1]=m2_[1][0]=sqrt(s3); m2_[0][2]=m2_[2][0]=sqrt(s2); m2_[1][2]=m2_[2][1]=sqrt(s1); // now compute the matrix element // amplitide complex amp = amplitude(-1); // epsilon piece Energy6 kin = (pow<4,1>(m1)*(-2*(sqr(m2) + sqr(m3)) + s1) + pow<4,1>(m2)*(-2*sqr(m3) + s2) + s3*(pow<4,1>(m3) + s1*s2 - sqr(m3)*(s1 + s2 + s3)) - sqr(m1)*(2*pow<4,1>(m2) + 2*pow<4,1>(m3) + sqr(m2)*(4*sqr(m3) - 3*(s1 + s2) - s3) + s1*(s1 + s2 + s3) - sqr(m3)*(3*s1 + s2 + 3*s3)) - sqr(m2)*(2*pow<4,1>(m3) + s2*(s1 + s2 + s3) - sqr(m3)*(s1 + 3*(s2 + s3))))/12.; return norm(amp*coupling_*GeV*GeV2)*kin/GeV2/GeV2/GeV2; } WidthCalculatorBasePtr VectorTo3PseudoScalarDalitz::threeBodyMEIntegrator(const DecayMode & ) const { int imode=0; // construct the integrator vector inweights; vector intype; vector inmass,inwidth; inweights.reserve(resonances().size()); intype.reserve(resonances().size()); inmass.reserve(resonances().size()); inwidth.reserve(resonances().size()); int iloc=-1; vector inpow(2,0.0); for(unsigned int ix=0;ixid); if(resonance) { ++iloc; inweights.push_back(weights()[iloc]); inmass.push_back(resonances()[ix]->mass); inwidth.push_back(abs(resonances()[ix]->width)); intype.push_back(resonances()[ix]->spectator+1); } } tcDecayIntegratorPtr decayer(this); WidthCalculatorBasePtr output( new_ptr(ThreeBodyAllOnCalculator (inweights,intype,inmass,inwidth,inpow, *this,imode, mode(0)->outgoing()[0]->mass(), mode(0)->outgoing()[1]->mass(), mode(0)->outgoing()[2]->mass()))); return output; } void VectorTo3PseudoScalarDalitz::dataBaseOutput(ofstream & output, bool header) const { if(header){output << "update decayers set parameters=\"";} output << "newdef " << name() << ":Coupling " << coupling_*GeV << "\n"; // parameters for the DalitzBase base class DalitzBase::dataBaseOutput(output,false); if(header){output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;} } complex VectorTo3PseudoScalarDalitz::resAmp(unsigned int i) const { // can't have a scalar here on spin/parity grounds assert(resonances()[i]->type%10!=1); // shouldn't have E691 stuff either assert(resonances()[i]->type%10!=1); // amplitude - static const Complex ii = Complex(0.,1.); Complex output = resonances()[i]->amp; if (resonances()[i]->type==ResonanceType::NonResonant) return output/GeV2; + // mass of the resonance + const Energy & mR = resonances()[i]->mass ; // locations of the outgoing particles const int &d1 = resonances()[i]->daughter1; const int &d2 = resonances()[i]->daughter2; const int &sp = resonances()[i]->spectator; // epsilon piece =eps(d1,d2,sp) double sign = (sp-d1)*(d1-d2)*(d2-sp)/2.; - // mass and width of the resonance - const Energy & mR = resonances()[i]->mass ; - const Energy & wR = resonances()[i]->width; - // momenta for the resonance decay - // off-shell - Energy pAB=sqrt(0.25*sqr(sqr(m2_[d1][d2]) -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/m2_[d1][d2]; - // on-shell - Energy pR=sqrt(0.25*sqr( mR*mR -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/mR; + // compute the Breit-Wigner times resonance form factor piece + output *= sign*resonances()[i]->BreitWigner(m2_[d1][d2],mOut_[d1],mOut_[d2]); // Blatt-Weisskopf factors - double fR=1, fD=1; - unsigned int power(1); // for the D decay Energy pD = sqrt(max(ZERO,(0.25*sqr(sqr(mD_)-sqr(mR)-sqr(mOut_[sp])) - sqr(mR*mOut_[sp]))/sqr(mD_))); Energy pDAB= sqrt( 0.25*sqr(sqr(mD_)-sqr(m2_[d1][d2])-sqr(mOut_[sp])) - sqr(m2_[d1][d2]*mOut_[sp]))/mD_; - double r1A(resonances()[i]->R*pR),r1B(resonances()[i]->R*pAB ); double r2A(parentRadius() *pD),r2B(parentRadius() *pDAB); // Blatt-Weisskopf factors and spin piece switch (resonances()[i]->type) { - case ResonanceType::Spin1: - fR=sqrt( (1. + sqr(r1A)) / (1. + sqr(r1B)) ); - fD=sqrt( (1. + sqr(r2A)) / (1. + sqr(r2B)) ); - power=3; - output *= fR*fD; + case ResonanceType::Spin1: case ResonanceType::Spin1GS : + output *= sqrt( (1. + sqr(r2A)) / (1. + sqr(r2B)) ); break; case ResonanceType::Spin2: - fR = sqrt( (9. + sqr(r1A)*(3.+sqr(r1A))) / (9. + sqr(r1B)*(3.+sqr(r1B)))); - fD = sqrt( (9. + sqr(r2A)*(3.+sqr(r2A))) / (9. + sqr(r2B)*(3.+sqr(r2B)))); - power=5; - output *= fR*fD; + output *= sqrt( (9. + sqr(r2A)*(3.+sqr(r2A))) / (9. + sqr(r2B)*(3.+sqr(r2B)))); // spin piece output *= (sqr(mD_) - sqr(mOut_[sp]) + sqr(m2_[d1][d2]))/(sqr(mD_)*sqr(m2_[d1][d2]))/GeV2* ((-sqr(m2_[sp][d1]) + sqr(m2_[sp][d2]))*sqr(m2_[d1][d2]) + (mD_ - mOut_[sp])*(mD_ + mOut_[sp])*(mOut_[d1] - mOut_[d2])*(mOut_[d1] + mOut_[d2])); break; default : assert(false); } - Energy gam = wR*pow(pAB/pR,power)*(mR/m2_[d1][d2])*fR*fR; - return sign*output/(sqr(mR)-sqr(m2_[d1][d2])-mR*gam*ii); + return output/GeV2; }