Page MenuHomeHEPForge
Files F7879083

No OneTemporary
Actions

View Options

diff --git a/Decay/ScalarMeson/CLEOD0toK0PipPim.cc b/Decay/ScalarMeson/CLEOD0toK0PipPim.cc
--- a/Decay/ScalarMeson/CLEOD0toK0PipPim.cc
+++ b/Decay/ScalarMeson/CLEOD0toK0PipPim.cc
@@ -1,131 +1,131 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the CLEOD0toK0PipPim class.
//
#include "CLEOD0toK0PipPim.h"
#include "ThePEG/Interface/ClassDocumentation.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"
using namespace Herwig;
CLEOD0toK0PipPim::CLEOD0toK0PipPim() : WeakDalitzDecay(5./GeV,1.5/GeV,true)
{}
IBPtr CLEOD0toK0PipPim::clone() const {
return new_ptr(*this);
}
IBPtr CLEOD0toK0PipPim::fullclone() const {
return new_ptr(*this);
}
void CLEOD0toK0PipPim::persistentOutput(PersistentOStream & ) const {
}
void CLEOD0toK0PipPim::persistentInput(PersistentIStream & , int) {
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<CLEOD0toK0PipPim,WeakDalitzDecay>
describeHerwigCLEOD0toK0PipPim("Herwig::CLEOD0toK0PipPim", "CLEOD0toK0PipPim.so");
void CLEOD0toK0PipPim::Init() {
static ClassDocumentation<CLEOD0toK0PipPim> documentation
("The CLEOD0toK0PipPim class implements the model of CLEO for"
" D0 -> Kbar0 pi+pi-",
"The CLEO fit of \\cite{Muramatsu:2002jp} was used for the decay $D^0\\to\\bar{K}^0\\pi^+\\pi^-$",
"\\bibitem{Muramatsu:2002jp} H.~Muramatsu {\\it et al.} "
"[CLEO Collaboration],Phys.\\ Rev.\\ Lett.\\ {\\bf 89} (2002) 251802"
"[Erratum-ibid.\\ {\\bf 90} (2003) 059901] [arXiv:hep-ex/0207067].\n");
}
void CLEOD0toK0PipPim::doinit() {
WeakDalitzDecay::doinit();
static const double degtorad = Constants::pi/180.;
// create the resonances
addResonance(DalitzResonance(getParticleData( 323),0.89166*GeV,0.0508*GeV,0,1,2,-0.11 , 321*degtorad));
addResonance(DalitzResonance(getParticleData( 113),0.7693 *GeV,0.1502*GeV,1,2,0,-1. , 0 ));
addResonance(DalitzResonance(getParticleData( 223),0.78257*GeV, 8.44*MeV,1,2,0,-0.037 , 114*degtorad));
addResonance(DalitzResonance(getParticleData( -323),0.89166*GeV,0.0508*GeV,0,2,1,-1.56 , 150*degtorad));
addResonance(DalitzResonance(getParticleData(9010221),0.977 *GeV,0.050 *GeV,1,2,0, 0.34 , 188*degtorad));
addResonance(DalitzResonance(getParticleData( 225),1.2754 *GeV,0.1851*GeV,1,2,0, 0.7 , 308*degtorad));
addResonance(DalitzResonance(getParticleData( 10211),1.310 *GeV,0.2720*GeV,1,2,0, 1.8 , 85*degtorad));
addResonance(DalitzResonance(getParticleData( -10321),1.412 *GeV,0.294 *GeV,0,2,1, 2.0 , 3*degtorad));
addResonance(DalitzResonance(getParticleData( -325),1.4256 *GeV,0.0985*GeV,0,2,1, 1.0 , 335*degtorad));
addResonance(DalitzResonance(getParticleData( -30323),1.717 *GeV,0.322 *GeV,0,2,1,-5.6 , 174*degtorad));
// D0 -> K- pi+ pi0
createMode(getParticleData(ParticleID::D0),
{getParticleData(ParticleID::Kbar0),
getParticleData(ParticleID::piplus),
getParticleData(ParticleID::piminus)});
}
void CLEOD0toK0PipPim::doinitrun() {
WeakDalitzDecay::doinitrun();
}
int CLEOD0toK0PipPim::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
int id0(parent->id());
// incoming particle must be D0
if(abs(id0)!=ParticleID::D0) return -1;
cc = id0==ParticleID::Dbar0;
// must be three decay products
if(children.size()!=3) return -1;
tPDVector::const_iterator pit = children.begin();
unsigned int npip(0),npim(0),nkm(0),nk0(0),npi0(0);
for( ;pit!=children.end();++pit) {
id0=(**pit).id();
if(id0 ==ParticleID::piplus) ++npip;
else if(id0 ==ParticleID::pi0) ++npi0;
else if(id0 ==ParticleID::piminus) ++npim;
else if(abs(id0)==ParticleID::K0) ++nk0;
else if(id0 ==ParticleID::K_L0) ++nk0;
else if(id0 ==ParticleID::K_S0) ++nk0;
else if(abs(id0)==ParticleID::Kplus) ++nkm;
}
if(npim==1&&npip==1&&nk0==1) return 0;
else return -1;
}
Complex CLEOD0toK0PipPim::amplitude(int ichan) const {
Complex output(0.);
static const Complex ii(0.,1.);
unsigned int imin=0, imax(resonances().size());
if(ichan>=0) {
imin=ichan;
imax=imin+1;
}
for(unsigned int ix=imin;ix<imax;++ix) {
// all resonances bar f0(980)
if(resonances()[ix].resonance->id()!=9010221) {
output += resAmp(ix);
}
// special treatment (Flatte) for f0(980)
else {
// output += resAmp(ix);
static const double gpi=0.09, gK=0.02;
const Energy & mAB = mInv(resonances()[ix].daughter1,resonances()[ix].daughter2);
Energy Gamma_pi = gpi*sqrt(0.25*sqr(mAB)-sqr(mOut(resonances()[ix].daughter1)));
Energy2 arg = 0.25*sqr(mAB)-sqr(mOut(resonances()[ix].spectator));
complex<Energy> Gamma_K = arg>=ZERO ? gK*sqrt(arg) : gK*ii*sqrt(-arg);
- output += resonances()[ix].amplitude*Complex(cos(resonances()[ix].phase),sin(resonances()[ix].phase))*GeV2/
+ output += resonances()[ix].amp*GeV2/
(sqr(resonances()[ix].mass)-sqr(mAB)-ii*resonances()[ix].mass*(Gamma_pi+Gamma_K));
}
}
if(ichan<0) output += 1.1*Complex(cos(5.93411945678072),sin(5.93411945678072));
return output;
}
diff --git a/Decay/ScalarMeson/WeakDalitzDecay.cc b/Decay/ScalarMeson/WeakDalitzDecay.cc
--- a/Decay/ScalarMeson/WeakDalitzDecay.cc
+++ b/Decay/ScalarMeson/WeakDalitzDecay.cc
@@ -1,159 +1,159 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the WeakDalitzDecay class.
//
#include "WeakDalitzDecay.h"
#include "ThePEG/Interface/ClassDocumentation.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"
using namespace Herwig;
WeakDalitzDecay::WeakDalitzDecay(InvEnergy rP, InvEnergy rR, bool useResonanceMass)
: rParent_(rP), rResonance_(rR), useResonanceMass_(useResonanceMass), maxWgt_(1.)
{}
void WeakDalitzDecay::persistentOutput(PersistentOStream & os) const {
os << resonances_ << maxWgt_ << weights_;
}
void WeakDalitzDecay::persistentInput(PersistentIStream & is, int) {
is >> resonances_ >> maxWgt_ >> weights_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractClass<WeakDalitzDecay,DecayIntegrator>
describeHerwigWeakDalitzDecay("Herwig::WeakDalitzDecay", "HwSMDecay.so");
void WeakDalitzDecay::Init() {
static ClassDocumentation<WeakDalitzDecay> documentation
("The WeakDalitzDecay class provides a base class for "
"weak three-body decays of bottom and charm mesons");
}
void WeakDalitzDecay::
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);
}
void WeakDalitzDecay::doinit() {
DecayIntegrator::doinit();
}
void WeakDalitzDecay::doinitrun() {
DecayIntegrator::doinitrun();
weights_.resize(mode(0)->channels().size());
maxWgt_ = mode(0)->maxWeight();
for(unsigned int iz=0;iz<mode(0)->channels().size();++iz) {
weights_[iz]=mode(0)->channels()[iz].weight();
}
}
double WeakDalitzDecay::me2(const int ichan, const Particle & part,
const tPDVector & ,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
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];
}
}
// compute the amplitude (in inherited class)
Complex amp = amplitude(ichan);
// now compute the matrix element
(*ME())(0,0,0,0)=amp;
return norm(amp);
}
void WeakDalitzDecay::createMode(tPDPtr in, tPDVector out) {
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,maxWgt_));
if(weights_.size()!=resonances_.size()) {
weights_=vector<double>(resonances_.size(),1./double(resonances_.size()));
}
unsigned int ix=0;
for(DalitzResonance res : resonances_) {
mode->addChannel((PhaseSpaceChannel(mode),0,res.resonance,0,res.spectator+1,1,res.daughter1+1,1,res.daughter2+1));
resetIntermediate(res.resonance,res.mass,res.width);
++ix;
}
addMode(mode);
}
Complex WeakDalitzDecay::resAmp(unsigned int i, bool gauss) const {
- Complex output = resonances_[i].amplitude*Complex(cos(resonances_[i].phase),sin(resonances_[i].phase));
+ Complex output = resonances_[i].amp;
// 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];
// on-shell
Energy pR=sqrt(0.25*sqr( mR*mR -sqr(mOut_[d1])-sqr(mOut_[d2])) - sqr(mOut_[d1]*mOut_[d2]))/mR;
// 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_;
// Blatt-Weisskopf factors
double fR=1, fD=1;
unsigned int power(1);
double r1A(rResonance_*pR),r1B(rResonance_*pAB );
double r2A(rParent_ *pD),r2B(rParent_ *pDAB);
// mass for thre denominator
Energy mDen = useResonanceMass_ ? mR : m2_[d1][d2];
// Blatt-Weisskopf factors and spin piece
switch (resonances_[i].resonance->iSpin()) {
case PDT::Spin0:
if(gauss) {
fR = exp(-(r1B-r1A)/12.);
fD = exp(-(r2B-r2A)/12.);
}
break;
case PDT::Spin1:
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) )/GeV2;
break;
case PDT::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/GeV2/GeV2*( 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
Energy gam = wR*pow(pAB/pR,power)*(mR/m2_[d1][d2])*fR*fR;
return output*GeV2/(sqr(mR)-sqr(m2_[d1][d2])-mR*gam*Complex(0.,1.));
}
diff --git a/Decay/ScalarMeson/WeakDalitzDecay.h b/Decay/ScalarMeson/WeakDalitzDecay.h
--- a/Decay/ScalarMeson/WeakDalitzDecay.h
+++ b/Decay/ScalarMeson/WeakDalitzDecay.h
@@ -1,304 +1,299 @@
// -*- C++ -*-
#ifndef Herwig_WeakDalitzDecay_H
#define Herwig_WeakDalitzDecay_H
//
// This is the declaration of the WeakDalitzDecay class.
//
#include "Herwig/Decay/DecayIntegrator.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
namespace Herwig {
using namespace ThePEG;
/**
* Struct to contain the properties of the intermediate
*/
struct DalitzResonance {
/**
* Default constructor
*/
DalitzResonance() {}
/**
* Constructor specifiying the parameters
*/
DalitzResonance(tPDPtr res, Energy m, Energy w,
unsigned int d1, unsigned int d2, unsigned int s,
- double amp, double phi)
+ double mag, double phi)
: resonance(res), mass(m),width(w),
daughter1(d1),daughter2(d2),spectator(s),
- amplitude(amp),phase(phi)
+ amp(mag*exp(Complex(0,phi)))
{}
/**
* Resonant particle
*/
tPDPtr resonance;
/**
* Mass of the resonance
*/
Energy mass;
/**
* Width of the resonance
*/
Energy width;
/**
* The children
*/
unsigned int daughter1;
unsigned int daughter2;
/**
* The spectactor
*/
unsigned int spectator;
/**
* The amplitude
*/
- double amplitude;
-
- /**
- * The phase
- */
- double phase;
+ Complex amp;
};
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
* @param x The resonance
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const DalitzResonance & x) {
os << x.resonance << ounit(x.mass,GeV) << ounit(x.width,GeV)
<< x.daughter1 << x.daughter2 << x.spectator
- << x.amplitude << x.phase;
+ << x.amp;
return os;
}
/**
* Input operator to allow the structure to be persistently read
* @param is The input stream
* @param x The resonance
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
DalitzResonance & x) {
is >> x.resonance >> iunit(x.mass,GeV) >> iunit(x.width,GeV)
>> x.daughter1 >> x.daughter2 >> x.spectator
- >> x.amplitude >> x.phase;
+ >> x.amp;
return is;
}
/**
* The WeakDalitzDecay class provides a base class for the implementation
* of weak three-body decays of bottom and charm mesons
*
* @see \ref WeakDalitzDecayInterfaces "The interfaces"
* defined for WeakDalitzDecay.
*/
class WeakDalitzDecay: public DecayIntegrator {
public:
/**
* The default constructor.
*/
WeakDalitzDecay(InvEnergy rP=5./GeV, InvEnergy rR=1.5/GeV, bool useResonanceMass=false);
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param outgoing The particles produced in the decay
* @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
double me2(const int ichan,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const;
/**
* Construct the SpinInfos for the particles produced in the decay
*/
virtual void constructSpinInfo(const Particle & part,
ParticleVector outgoing) 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);
//@}
/**
* 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:
/**
* Add a new resonance
*/
void addResonance(const DalitzResonance & R) {resonances_.push_back(R);}
/**
* Set up the phase-space
*/
void createMode(tPDPtr in, tPDVector out);
/**
* Calculate the amplitude
*/
virtual Complex amplitude(int ichan) const = 0;
/**
* Calculate the amplitude for the ith resonance
*/
Complex resAmp(unsigned int i,bool gauss=false) const;
/**
* Access to the resonances
*/
const vector<DalitzResonance> & resonances() const {return resonances_;}
/**
* Access to the invariants
*/
const Energy & mInv(unsigned int i,unsigned int j) const {
return m2_[i][j];
}
/**
* Masses of the children
*/
const Energy & mOut(unsigned int i) const {
return mOut_[i];
}
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();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
WeakDalitzDecay & operator=(const WeakDalitzDecay &);
private:
/**
* The radii for the Blatt-Weisskopf form factors
*/
//@{
/**
* For the decaying particles
*/
InvEnergy rParent_;
/**
* For the intermediate resonances
*/
InvEnergy rResonance_;
//@}
/**
* Vector containing the intermediate resonances
*/
vector<DalitzResonance> resonances_;
/**
* Choice of the mass to use in the denominator of the expressions
*/
bool useResonanceMass_;
private:
/**
* Parameters for the phase-space sampling
*/
//@{
/**
* Maximum weight for the decay
*/
double maxWgt_;
/**
* Weights for the phase-space channels
*/
vector<double> weights_;
//@}
private :
/**
* Storage of the kinematics
*/
//@{
/**
* Mass of the parent
*/
mutable Energy mD_;
/**
* Masses of the children
*/
mutable Energy mOut_[3];
/**
* Masses of the children
*/
mutable Energy m2_[3][3];
//@}
/**
* Spin density matrix
*/
mutable RhoDMatrix _rho;
};
}
#endif /* Herwig_WeakDalitzDecay_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:30 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805834
Default Alt Text
(18 KB)

Event Timeline