Page MenuHomeHEPForge
Files F8308762

No OneTemporary
Actions

View Options

diff --git a/Decay/Dalitz/ScalarTo3ScalarDalitz.cc b/Decay/Dalitz/ScalarTo3ScalarDalitz.cc
--- a/Decay/Dalitz/ScalarTo3ScalarDalitz.cc
+++ b/Decay/Dalitz/ScalarTo3ScalarDalitz.cc
@@ -1,215 +1,215 @@
// -*- 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;
IBPtr ScalarTo3ScalarDalitz::clone() const {
return new_ptr(*this);
}
IBPtr ScalarTo3ScalarDalitz::fullclone() const {
return new_ptr(*this);
}
void ScalarTo3ScalarDalitz::persistentOutput(PersistentOStream & os) const {
os << f0gpi_ << f0gK_ << useResonanceMass_ ;
}
void ScalarTo3ScalarDalitz::persistentInput(PersistentIStream & is, int) {
is >> f0gpi_ >> f0gK_ >> useResonanceMass_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ScalarTo3ScalarDalitz,DalitzBase>
describeHerwigScalarTo3ScalarDalitz("Herwig::ScalarTo3ScalarDalitz", "HwDalitzDecay.so");
void ScalarTo3ScalarDalitz::Init() {
static ClassDocumentation<ScalarTo3ScalarDalitz> documentation
("The ScalarTo3ScalarDalitz class provides a base class for "
"weak three-body decays of bottom and charm mesons");
static Parameter<ScalarTo3ScalarDalitz,double> interfacegPi
("f0gPi",
"The g_pi coupling for the f_0(980) width",
- &ScalarTo3ScalarDalitz::f0gpi_, 0.09, 0.0, 1.,
+ &ScalarTo3ScalarDalitz::f0gpi_, 0.09, 0.0, 10.,
false, false, Interface::limited);
static Parameter<ScalarTo3ScalarDalitz,double> interfacegK
("f0gK",
"The g_K coupling for the f_0(980) width",
- &ScalarTo3ScalarDalitz::f0gK_, 0.02, 0.0, 1.,
+ &ScalarTo3ScalarDalitz::f0gK_, 0.02, 0.0, 10.,
false, false, Interface::limited);
static Switch<ScalarTo3ScalarDalitz,bool> 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<tPPtr>(&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<Lorentz5Momentum> & 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<tPPtr>(&part),incoming);
}
// set the kinematics
mD_ = part.mass();
for(unsigned int ix=0;ix<momenta.size();++ix) {
mOut_[ix]=momenta[ix].mass();
for(unsigned int iy=ix+1;iy<momenta.size();++iy) {
m2_[ix][iy]=(momenta[ix]+momenta[iy]).m();
m2_[iy][ix]=m2_[ix][iy];
}
}
// now compute the matrix element
Complex amp = amplitude(ichan);
(*ME())(0,0,0,0) = amp;
return norm(amp);
}
Complex ScalarTo3ScalarDalitz::resAmp(unsigned int i) const {
static const Complex ii = Complex(0.,1.);
Complex output = resonances()[i].amp;
if (resonances()[i].type==ResonanceType::NonResonant) return output;
// 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;
// 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];
if(resonances()[i].type==ResonanceType::BABARf0) {
double rho = 2.*pAB/m2_[d1][d2];
output *= GeV2/(sqr(resonances()[i].mass)-sqr(m2_[d1][d2])-
ii*resonances()[i].mass*resonances()[i].width*rho);
return output;
}
else if (resonances()[i].type==ResonanceType::Flattef0) {
- Energy mpi = getParticleData(211)->mass();
+ Energy mpi = getParticleData(111)->mass();
Energy mK = getParticleData(321)->mass();
Energy Gamma_pi = f0gpi_*sqrt(0.25*sqr(m2_[d1][d2])-sqr(mpi));
Energy2 arg = 0.25*sqr(m2_[d1][d2])-sqr(mK);
complex<Energy> Gamma_K = arg>=ZERO ? f0gK_*sqrt(arg) : f0gK_*ii*sqrt(-arg);
output *= GeV2/(sqr(resonances()[i].mass)-sqr(m2_[d1][d2])-ii*resonances()[i].mass*(Gamma_pi+Gamma_K));
return output;
}
else if (resonances()[i].type==ResonanceType::Spin0Complex) {
complex<Energy> sR(resonances()[i].mass,resonances()[i].width);
output *= GeV2/(sqr(sR)-sqr(m2_[d1][d2]));
return output;
}
// on-shell
Energy pR=sqrt(0.25*sqr( mR*mR -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/mR;
// Blatt-Weisskopf factors
double fR=1, fD=1;
unsigned int power(1);
if(resonances()[i].type!=ResonanceType::Spin0 &&
resonances()[i].type!=ResonanceType::Spin0E691) {
// 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);
// mass for thre denominator
Energy mDen = useResonanceMass_ ? mR : 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);
}
// Blatt-Weisskopf factors and spin piece
switch (resonances()[i].type) {
case ResonanceType::Spin0Gauss:
fR = exp(-(r1B-r1A)/12.);
fD = exp(-(r2B-r2A)/12.);
break;
case ResonanceType::Spin1: case ResonanceType::Spin1E691 :
fR=sqrt( (1. + sqr(r1A)) / (1. + sqr(r1B)) );
fD=sqrt( (1. + sqr(r2A)) / (1. + sqr(r2B)) );
power=3;
output *= fR*fD*(sqr(m2_[d2][sp])-sqr(m2_[d1][sp])
+ ( sqr(mD_)-sqr(mOut_[sp]))*(sqr(mOut_[d1])-sqr(mOut_[d2]))/sqr(mDen) )/denom;
break;
case ResonanceType::Spin2: case ResonanceType::Spin2E691:
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/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.);
break;
default :
assert(false);
}
}
// multiply by Breit-Wigner piece and return
if (resonances()[i].type/10 == 1 ) {
return output*sqrt(0.5*wR/GeV/Constants::pi)*GeV/(m2_[d1][d2]-mR-complex<Energy>(ZERO,0.5*wR));
}
else {
Energy gam = wR*pow(pAB/pR,power)*(mR/m2_[d1][d2])*fR*fR;
return output*GeV2/(sqr(mR)-sqr(m2_[d1][d2])-mR*gam*ii);
}
}
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 << "newdef " << name() << ":f0gPi " << f0gpi_ << "\n";
output << "newdef " << name() << ":f0gK " << f0gK_ << "\n";
output << "\n";
if(header) {
output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:05 PM (1 d, 42 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4022885
Default Alt Text
(9 KB)

Event Timeline