Page MenuHomeHEPForge
Files F7877142

No OneTemporary
Actions

View Options

diff --git a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc
--- a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc
+++ b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.cc
@@ -1,362 +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;ix<incoming_.size();++ix) {
if(mode(ix)) maxweight_[ix] = mode(ix)->maxWeight();
}
}
}
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;ix<incoming_.size();++ix) {
tPDPtr in = getParticleData(incoming_[ix]);
tPDVector out = {getParticleData(outgoingf_[ix]),
getParticleData(outgoinga_[ix])};
if(in&&out[0]&&out[1])
mode = new_ptr(PhaseSpaceMode(in,out,maxweight_[ix]));
else
mode=PhaseSpaceModePtr();
addMode(mode);
}
}
VectorMeson2BaryonsDecayer::VectorMeson2BaryonsDecayer()
: gm_ ({0.00163222616377,0.00158881446341,0.00163819673075,0.000944247071286,0.00141162244048 ,0.00150424724997 ,0.00093350341391,0.00107528142563,0.000860759359361,0.00103222902272,0.00098818310509,0.00119542938517 ,0.000941856299908,0.00108249234586 ,0.00102912698738 ,0.000561082932891,0.000526423800011,0.000500391020504,0.000678994938986}),
ge_ ({0.00135736265293,0.0015117595557 ,0.00136696938558,0.001988067505 ,0.000852659560453,0.000801719253474,0.00150640470181,0.00153402258271,0.0015323974485 ,0. ,0.00084599445285,0.000621041230885,0.000599393812308,0.000327664954147,0.000664382626732,0.000260326162249,0.000362882855521,0.00025226758188 ,0.000397805035356}),
phi_ ({0. ,0. ,0.740 ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. }),
incoming_ ({443 ,443 ,443 ,443 ,443 ,443 ,443 ,443 ,443 ,100443 ,100443 ,100443 ,100443 ,100443 ,100443 ,100443 ,100443 ,100443 ,100443 }),
outgoingf_({ 2212 , 2112 , 3122 , 3212 , 3312 , 3322 , 3114 , 3224 , 3214 , 2212 , 2112 , 3122 , 3212 , 3312 , 3322 , 3114 , 3224 , 3214 , 3314 }),
outgoinga_({-2212 ,-2112 ,-3122 ,-3212 ,-3312 ,-3322 ,-3114 ,-3224 ,-3214 ,-2212 ,-2112 ,-3122 ,-3212 ,-3312 ,-3322 ,-3114 ,-3224 ,-3214 ,-3314 }),
maxweight_({1.6 ,1.6 ,2.1 ,1.5 ,2. ,2.1 ,7. ,10. ,6.5 ,1.7 ,1.7 ,2.5 ,2.5 ,2. ,2. ,30. ,35. ,41. , 0.001 })
{}
int VectorMeson2BaryonsDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
if(children.size()!=2) return -1;
int id(parent->id());
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<incoming_.size());
return imode;
}
IBPtr VectorMeson2BaryonsDecayer::clone() const {
return new_ptr(*this);
}
IBPtr VectorMeson2BaryonsDecayer::fullclone() const {
return new_ptr(*this);
}
void VectorMeson2BaryonsDecayer::persistentOutput(PersistentOStream & os) const {
os << gm_ << ge_ << phi_ << incoming_ << outgoingf_ << outgoinga_ << maxweight_;
}
void VectorMeson2BaryonsDecayer::persistentInput(PersistentIStream & is, int) {
is >> gm_ >> ge_ >> phi_ >> incoming_ >> outgoingf_ >> outgoinga_ >> maxweight_;
}
DescribeClass<VectorMeson2BaryonsDecayer,DecayIntegrator>
describeHerwigVectorMeson2BaryonsDecayer("Herwig::VectorMeson2BaryonsDecayer", "HwVMDecay.so");
void VectorMeson2BaryonsDecayer::Init() {
static ClassDocumentation<VectorMeson2BaryonsDecayer> documentation
("The VectorMeson2BaryonsDecayer class is designed for "
"the decay of vector mesons to baryons.");
static ParVector<VectorMeson2BaryonsDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&VectorMeson2BaryonsDecayer::incoming_,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMeson2BaryonsDecayer,int> interfaceOutcoming1
("OutgoingBaryons",
"The PDG code for the outgoing fermion",
&VectorMeson2BaryonsDecayer::outgoingf_,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMeson2BaryonsDecayer,int> interfaceOutcoming2
("OutgoingAntiBaryons",
"The PDG code for the second outgoing anti-fermion",
&VectorMeson2BaryonsDecayer::outgoinga_,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMeson2BaryonsDecayer,double> interfaceGM
("GM",
"The value of the GM form factor",
&VectorMeson2BaryonsDecayer::gm_, -1, 0., -1000., 1000.,
false, false, Interface::limited);
static ParVector<VectorMeson2BaryonsDecayer,double> interfaceGE
("GE",
"The value of the GE form factor",
&VectorMeson2BaryonsDecayer::ge_, -1, 0., -1000., 1000.,
false, false, Interface::limited);
static ParVector<VectorMeson2BaryonsDecayer,double> interfacePhi
("Phi",
"The phase of the GE form factor",
&VectorMeson2BaryonsDecayer::phi_, -1, 0., -Constants::pi, Constants::pi,
false, false, Interface::limited);
static ParVector<VectorMeson2BaryonsDecayer,double> 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<tPPtr>(&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<Lorentz5Momentum> & 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<tPPtr>(&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
- //double mesum(0.);
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);
- //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 << " " << 1./3.*(40.*norm(GM)/9.+16.*sqr(outgoing[0]->mass()/part.mass())*norm(GE)) << "\n";
+ // 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;
+ 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<incoming_.size();++ix) {
if(ix<initsize_) {
output << "newdef " << name() << ":Incoming " << ix << " "
<< incoming_[ix] << "\n";
output << "newdef " << name() << ":OutgoingFermion " << ix << " "
<< outgoingf_[ix] << "\n";
output << "newdef " << name() << ":OutgoingAntiFermion " << ix << " "
<< outgoinga_[ix] << "\n";
output << "newdef " << name() << ":GM " << ix << " "
<< gm_[ix] << "\n";
output << "newdef " << name() << ":GE " << ix << " "
<< ge_[ix] << "\n";
output << "newdef " << name() << ":Phi " << ix << " "
<< phi_[ix] << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< maxweight_[ix] << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< incoming_[ix] << "\n";
output << "insert " << name() << ":OutgoingFermion " << ix << " "
<< outgoingf_[ix] << "\n";
output << "insert " << name() << ":OutgoingAntiFermion " << ix << " "
<< outgoinga_[ix] << "\n";
output << "insert " << name() << ":GM " << ix << " "
<< gm_[ix] << "\n";
output << "insert " << name() << ":GE " << ix << " "
<< ge_[ix] << "\n";
output << "insert " << name() << ":Phi " << ix << " "
<< phi_[ix] << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< maxweight_[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/PDT/BaryonWidthGenerator.cc b/PDT/BaryonWidthGenerator.cc
--- a/PDT/BaryonWidthGenerator.cc
+++ b/PDT/BaryonWidthGenerator.cc
@@ -1,259 +1,256 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the BaryonWidthGenerator class.
//
#include "BaryonWidthGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/Baryon/Baryon1MesonDecayerBase.h"
using namespace Herwig;
using namespace ThePEG;
void BaryonWidthGenerator::persistentOutput(PersistentOStream & os) const {
os << _baryondecayers << _modeloc;
}
void BaryonWidthGenerator::persistentInput(PersistentIStream & is, int) {
is >> _baryondecayers >> _modeloc;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<BaryonWidthGenerator,GenericWidthGenerator>
describeHerwigBaryonWidthGenerator("Herwig::BaryonWidthGenerator", "HwBaryonDecay.so");
void BaryonWidthGenerator::Init() {
static ClassDocumentation<BaryonWidthGenerator> documentation
("The BaryonWidthGenerator class is designed for the calculation of the"
" running width for baryons.");
static RefVector<BaryonWidthGenerator,Baryon1MesonDecayerBase> interfaceBaryonDecayers
("BaryonDecayers",
"Pointers to the baryon decayers to get the couplings",
&BaryonWidthGenerator::_baryondecayers, -1, false, false, true, true, false);
static ParVector<BaryonWidthGenerator,int> interfaceModeLocation
("ModeLocation",
"The location of the mode in the decayer",
&BaryonWidthGenerator::_modeloc, 0, -1, 0, 0,
false, false, false);
}
void BaryonWidthGenerator::setupMode(tcDMPtr mode, tDecayIntegratorPtr decayer,
unsigned int) {
// cast the decayer
tBaryon1MesonDecayerBasePtr
baryon(dynamic_ptr_cast<tBaryon1MesonDecayerBasePtr>(decayer));
if(baryon) {
int dmode(baryon->findMode(*mode));
if(dmode<0) {
_baryondecayers.push_back(Baryon1MesonDecayerBasePtr());
_modeloc.push_back(-1);
return;
}
else {
_baryondecayers.push_back(baryon);
_modeloc.push_back(dmode);
}
}
else {
_baryondecayers.push_back(Baryon1MesonDecayerBasePtr());
_modeloc.push_back(-1);
}
}
void BaryonWidthGenerator::dataBaseOutput(ofstream & output, bool header) {
if(header) output << "update Width_Generators set parameters=\"";
// info from the base class
GenericWidthGenerator::dataBaseOutput(output,false);
// info from this class
for(unsigned int ix=0;ix<_baryondecayers.size();++ix) {
if(_baryondecayers[ix]) {
output << "insert " << name() << ":BaryonDecayers " << ix
<< " " << _baryondecayers[ix]->fullName() << "\n";
}
else {
output << "insert " << name() << ":BaryonDecayers " << ix
<< " NULL \n";
}
output << "insert " << name() << ":ModeLocation " << ix
<< " " << _modeloc[ix] << "\n";
}
if(header) output << "\n\" where BINARY ThePEGName=\"" << name() << "\";"
<< endl;
}
Energy BaryonWidthGenerator::partial2BodyWidth(int imode, Energy q,Energy m1,
Energy m2) const {
using Constants::pi;
Complex A1,A2,A3,B1,B2,B3;
if(q<m1+m2)
return ZERO;
// mode from the base class
int mecode(MEcode(imode));
if(mecode<=100){
return GenericWidthGenerator::partial2BodyWidth(imode,q,m1,m2);
}
// calcluate the decay momentum
Energy2 q2(q*q),m12(m1*m1),m22(m2*m2),
pcm2(0.25*(q2*(q2-2.*m12-2.*m22)+(m12-m22)*(m12-m22))/q2);
Energy pcm(sqrt(pcm2)),gam(ZERO),msum(q+m1);
// Energy m0(mass());
Energy2 fact1((q+m1)*(q+m1)-m22),fact2((q-m1)*(q-m1)-m22),fact3(q2+m12-m22);
// 1/2 -> 1/2 0
if(mecode==101) {
_baryondecayers[imode]->halfHalfScalarCoupling(_modeloc[imode],q,m1,m2,A1,B1);
double Afact1((A1*conj(A1)).real()),Bfact1((B1*conj(B1)).real());
gam = 0.125/pi/q2*pcm*(Afact1*fact1+Bfact1*fact2);
/*
Energy Qp(sqrt(pow(q+m1,2)-pow(m2,2))),Qm(sqrt(pow(q-m1,2)-pow(m2,2)));
Complex h1(2.*Qp*A1),h2(-2.*Qm*B1);
cout << "testing 1/2->1/2 0 "
<< gam << " "
<< real(h1*conj(h1)+h2*conj(h2))/32./pi*pcm/q2 << " "
<< real(h1*conj(h1)+h2*conj(h2))/32./pi*pcm/q2/gam << endl;
*/
}
// 1/2 -> 1/2 1
else if(mecode==102) {
_baryondecayers[imode]->halfHalfVectorCoupling(_modeloc[imode],q,m1,m2,
A1,A2,B1,B2);
double
Afact1((A1*conj(A1)).real()),Afact2((A2*conj(A2)).real()),
Bfact1((B1*conj(B1)).real()),Bfact2((B2*conj(B2)).real()),
Afact4((A1*conj(A2)+A2*conj(A1)).real()),
Bfact4((B1*conj(B2)+B2*conj(B1)).real());
Energy2 me(2.*(fact1*Bfact1+fact2*Afact1));
if(m2>1e-10*GeV) {
me+=1./m22*(fact1*Bfact1*(q-m1)*(q-m1)-2.*q*pcm*Bfact4*q*pcm*(q-m1)/msum
+fact2*q2*Bfact2*pcm2/msum/msum
+fact2*Afact1*msum *msum +2.*q*pcm*Afact4*q*pcm
+fact1*q2*Afact2*pcm2/msum/msum);
}
gam = pcm/8./pi/q2*me;
// test of the matrix element
/*
Energy2 Qp(sqrt(pow(q+m1,2)-pow(m2,2))),Qm(sqrt(pow(q-m1,2)-pow(m2,2)));
double r2(sqrt(2.));
Complex h1(2.*r2*Qp*B1),h2(-2.*r2*Qm*A1),h3(0.),h4(0.);
if(m2>1e-10*GeV)
{
h3=2./m2*(Qp*(q-m1)*B1-Qm*q*B2*pcm/(q+m1));
h4=2./m2*(Qm*(q+m1)*A1+Qp*q*A2*pcm/(q+m1));
}
cout << "testing 1/2->1/2 0 "
<< gam << " "
<< real(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+h4*conj(h4))/32./pi*pcm/q2
<< " "
<< real(h1*conj(h1)+h2*conj(h2)+h3*conj(h3)+h4*conj(h4))/32./pi*pcm/q2/gam
<< endl;
*/
}
// 1/2 -> 3/2 0
else if(mecode==103) {
_baryondecayers[imode]->halfThreeHalfScalarCoupling(_modeloc[imode],q,m1,m2,A1,B1);
double Afact1((A1*conj(A1)).real()),Bfact1((B1*conj(B1)).real());
gam = 0.25/pi/3./msum/msum/m12*pcm*pcm2*(Afact1*fact1+Bfact1*fact2);
/*
Energy2 Qp(sqrt(pow(q+m1,2)-pow(m2,2))),Qm(sqrt(pow(q-m1,2)-pow(m2,2)));
double r23(sqrt(2./3.));
Complex h1(-2.*r23*pcm*q/m1*Qm*B1/(q+m1)),h2( 2.*r23*pcm*q/m1*Qp*A1/(q+m1));
cout << "testing 1/2->3/2 0 "
<< gam << " "
<< real(h1*conj(h1)+h2*conj(h2))/32./pi*pcm/q2 << " "
<< real(h1*conj(h1)+h2*conj(h2))/32./pi*pcm/q2/gam << endl;
*/
}
// 1/2 -> 3/2 1
else if(mecode==104) {
Energy Qp(sqrt(fact1)),Qm(sqrt(fact2));
double r2(sqrt(2.)),r3(sqrt(3.));
complex<Energy> h1(-2.*Qp*A1),h2(2.*Qm*B1),h5(ZERO),h6(ZERO);
complex<Energy> h3(-2./r3*Qp*(A1-fact2/m1*A2/msum));
complex<Energy> h4( 2./r3*Qm*(B1-fact1/m1*B2/msum));
if(m2>1e-10*GeV) {
h5=-2.*r2/r3/m1/m2*Qp*(0.5*(q2-m12-m22)*A1+0.5*fact2*(q+m1)*A2/msum
+q2*pcm*pcm*A3/sqr(msum));
h6= 2.*r2/r3/m1/m2*Qm*(0.5*(q2-m12-m22)*B1-0.5*fact1*(q-m1)*B2/msum
+q2*pcm*pcm*B3/sqr(msum));
}
gam=real(+h1*conj(h1)+h2*conj(h2)+h3*conj(h3)
+h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/32./pi*pcm/q2;
}
// 3/2 -> 1/2 0
else if(mecode==105) {
_baryondecayers[imode]->threeHalfHalfScalarCoupling(_modeloc[imode],q,m1,m2,A1,B1);
double Afact1((A1*conj(A1)).real()),Bfact1((B1*conj(B1)).real());
gam = 0.125/3./pi/msum/msum/q2*pcm*pcm2*(Afact1*fact1+Bfact1*fact2);
- /*
- Energy2 Qp(sqrt(pow(q+m1,2)-pow(m2,2))),Qm(sqrt(pow(q-m1,2)-pow(m2,2)));
- double r23(sqrt(2./3.));
- Complex h1(-2.*r23*pcm*Qm*B1/(q+m1)),
- h2( 2.*r23*pcm*Qp*A1/(q+m1));
- cout << "testing 3/2->1/2 0 "
- << gam << " "
- << real(h1*conj(h1)+h2*conj(h2))/64./pi*pcm/q2 << " "
- << real(h1*conj(h1)+h2*conj(h2))/64./pi*pcm/q2/gam << endl;
- */
+ // Energy Qp(sqrt(sqr(q+m1)-sqr(m2))),Qm(sqrt(sqr(q-m1)-sqr(m2)));
+ // double r23(sqrt(2./3.));
+ // complex<Energy> h1(-2.*r23*pcm*Qm*B1/(q+m1)), h2( 2.*r23*pcm*Qp*A1/(q+m1));
+ // generator()->log() << "testing 3/2->1/2 0 "
+ // << gam/GeV << " "
+ // << real(h1*conj(h1)+h2*conj(h2))/64./pi*pcm/q2/GeV << " "
+ // << real(h1*conj(h1)+h2*conj(h2))/64./pi*pcm/q2/gam << endl;
}
// 3/2 -> 1/2 1
else if(mecode==106) {
_baryondecayers[imode]->threeHalfHalfVectorCoupling(_modeloc[imode],q,m1,m2,
A1,A2,A3,B1,B2,B3);
Energy Qp(sqrt(fact1)),Qm(sqrt(fact2));
double r2(sqrt(2.)),r3(sqrt(3.));
complex<Energy> h1(-2.*Qp*A1),h2(2.*Qm*B1),h5(ZERO),h6(ZERO);
complex<Energy> h3(-2./r3*Qp*(A1-fact2/q*(A2/msum)));
complex<Energy> h4( 2./r3*Qm*(B1-fact1/q*(B2/msum)));
if(m2>1e-10*GeV) {
h5=-2.*r2/r3/q/m2*Qp*(0.5*(m12-q2-m22)*A1+0.5*fact2*(m1+q)*A2/msum
+q2*pcm*pcm*(A3/sqr(msum)));
h6= 2.*r2/r3/q/m2*Qm*(0.5*(m12-q2-m22)*B1-0.5*fact1*(m1-q)*(B2/msum)
+q2*pcm*pcm*(B3/sqr(msum)));
}
gam=real(+h1*conj(h1)+h2*conj(h2)+h3*conj(h3)
+h4*conj(h4)+h5*conj(h5)+h6*conj(h6))/64./pi*pcm/q2;
}
// 3/2 -> 3/2 0
else if(mecode==107) {
_baryondecayers[imode]->threeHalfThreeHalfScalarCoupling(_modeloc[imode],q,m1,m2,
A1,A2,B1,B2);
complex<InvEnergy2> A2byM2 = A2/sqr(msum);
complex<InvEnergy2> B2byM2 = B2/sqr(msum);
double Afact1((A1*conj(A1)).real());
InvEnergy4 Afact2((A2byM2*conj(A2byM2)).real());
double Bfact1((B1*conj(B1)).real());
InvEnergy4 Bfact2((B2byM2*conj(B2byM2)).real());
InvEnergy2 Afact4((A1*conj(A2byM2)+A2byM2*conj(A1)).real());
InvEnergy2 Bfact4((B1*conj(B2byM2)+B2byM2*conj(B1)).real());
gam = pcm/36./pi/q2/q2/m12*
( fact1*(Afact2*q2*q2*pcm2*pcm2
+0.25*Afact1*(fact3*fact1+10.*q2*m12)
+0.5*Afact4*q2*pcm*pcm*(fact3+q*m1))+
fact2*(Bfact2*q2*q2*pcm2*pcm2
+0.25*Bfact1*(fact3*fact2+10.*q2*m12)
+0.5*Bfact4*q2*pcm*pcm*(fact3-q*m1)));
}
else {
throw Exception() << "Unknown type of mode " << mecode
<< " in BaryonWidthGenerator::partial2BodyWidth() "
<< Exception::abortnow;
}
return gam*MEcoupling(imode)*MEcoupling(imode);
}
void BaryonWidthGenerator::doinit() {
if(initialize()) {
_baryondecayers.clear();
_modeloc.clear();
}
GenericWidthGenerator::doinit();
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 2:51 PM (1 d, 12 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3804841
Default Alt Text
(28 KB)

Event Timeline