Page MenuHomeHEPForge
Files F7879065

No OneTemporary
Actions

View Options

diff --git a/Decay/VectorMeson/VectorMesonVectorPScalarDecayer.cc b/Decay/VectorMeson/VectorMesonVectorPScalarDecayer.cc
--- a/Decay/VectorMeson/VectorMesonVectorPScalarDecayer.cc
+++ b/Decay/VectorMeson/VectorMesonVectorPScalarDecayer.cc
@@ -1,472 +1,474 @@
// -*- C++ -*-
//
// VectorMesonVectorPScalarDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VectorMesonVectorPScalarDecayer class.
//
#include "VectorMesonVectorPScalarDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Helicity/VectorSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void VectorMesonVectorPScalarDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<_incoming.size();++ix) {
if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight();
}
}
}
void VectorMesonVectorPScalarDecayer::doinit() {
DecayIntegrator::doinit();
// check consistency of the parameters
unsigned int isize=_incoming.size();
if(isize!=_outgoingV.size()||isize!=_outgoingP.size()||
isize!=_maxweight.size()||isize!=_coupling.size())
throw InitException() << "Inconsistent parameters in "
<< "VectorMesonVectorPScalarDecayer" << 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(_outgoingV[ix]),
getParticleData(_outgoingP[ix])};
if(in&&out[0]&&out[1])
mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix]));
else
mode=PhaseSpaceModePtr();
addMode(mode);
}
}
VectorMesonVectorPScalarDecayer::VectorMesonVectorPScalarDecayer()
: _coupling(73), _incoming(73), _outgoingV(73), _outgoingP(73), _maxweight(73) {
// intermediates
generateIntermediates(false);
// rho -> gamma pi modes
_incoming[0] = 113; _outgoingV[0] = 22; _outgoingP[0] = 111;
_coupling[0] = 0.2527/GeV; _maxweight[0] = 1.7;
_incoming[1] = 213; _outgoingV[1] = 22; _outgoingP[1] = 211;
_coupling[1] = 0.2210/GeV; _maxweight[1] = 1.7;
// rho -> gamma eta mode
_incoming[2] = 113; _outgoingV[2] = 22; _outgoingP[2] = 221;
_coupling[2] = 0.492/GeV; _maxweight[2] = 1.7;
// omega -> gamma pi
_incoming[3] = 223; _outgoingV[3] = 22; _outgoingP[3] = 111;
_coupling[3] = 0.7279465/GeV; _maxweight[3] = 1.7;
// omega -> gamma eta
_incoming[4] = 223; _outgoingV[4] = 22; _outgoingP[4] = 221;
_coupling[4] = 0.143/GeV; _maxweight[4] = 1.7;
// phi -> gamma pi
_incoming[5] = 333; _outgoingV[5] = 22; _outgoingP[5] = 111;
_coupling[5] = 0.0397/GeV; _maxweight[5] = 1.7;
// phi -> gamma eta
_incoming[6] = 333; _outgoingV[6] = 22; _outgoingP[6] = 221;
_coupling[6] = 0.212/GeV; _maxweight[6] = 1.7;
// phi -> gamma eta'
_incoming[7] = 333; _outgoingV[7] = 22; _outgoingP[7] = 331;
_coupling[7] = 0.219/GeV; _maxweight[7] = 1.8;
// phi -> omega pi
_incoming[8] = 333; _outgoingV[8] = 223; _outgoingP[8] = 111;
_coupling[8] = 0.0417/GeV; _maxweight[8] = 1.9;
// phi' -> K* K
_incoming[9] = 100333; _outgoingV[9] = 323; _outgoingP[9] = -321;
_coupling[9] = 3.934/GeV; _maxweight[9] = 5.;
_incoming[10] = 100333; _outgoingV[10] = 313; _outgoingP[10] = -311;
_coupling[10] = 4.011/GeV; _maxweight[10] = 5.;
// K* -> gamma K
_incoming[11] = 313; _outgoingV[11] = 22; _outgoingP[11] = 311;
_coupling[11] = 0.384/GeV; _maxweight[11] = 1.7;
_incoming[12] = 323; _outgoingV[12] = 22; _outgoingP[12] = 321;
_coupling[12] = 0.253/GeV; _maxweight[12] = 1.7;
// d* decay
_incoming[13] = 423; _outgoingV[13] = 22; _outgoingP[13] = 421;
_coupling[13] = 0.616/GeV; _maxweight[13] = 1.7;
_incoming[14] = 413; _outgoingV[14] = 22; _outgoingP[14] = 411;
_coupling[14] = 0.152/GeV; _maxweight[14] = 1.7;
// D_s* decays
_incoming[15] = 433; _outgoingV[15] = 22; _outgoingP[15] = 431;
_coupling[15] = 0.764/GeV; _maxweight[15] = 1.7;
// B_s* decays
_incoming[16] = 533; _outgoingV[16] = 22; _outgoingP[16] = 531;
_coupling[16] = 0.248/GeV; _maxweight[16] = 1.7;
// B_c* decays
_incoming[17] = 543; _outgoingV[17] = 22; _outgoingP[17] = 541;
_coupling[17] = 0.266/GeV; _maxweight[17] = 1.7;
// B* decay
_incoming[18] = 523; _outgoingV[18] = 22; _outgoingP[18] = 521;
_coupling[18] = 0.553/GeV; _maxweight[18] = 1.7;
_incoming[19] = 513; _outgoingV[19] = 22; _outgoingP[19] = 511;
_coupling[19] = 0.310/GeV; _maxweight[19] = 1.7;
// rho'' eta rho
_incoming[20] = 30113; _outgoingV[20] = 113; _outgoingP[20] = 221;
_coupling[20] = 2.663/GeV; _maxweight[20] = 4.;
_incoming[21] = 30213; _outgoingV[21] = 213; _outgoingP[21] = 221;
_coupling[21] = 2.663/GeV; _maxweight[21] = 4.;
// rho '' K* K
_incoming[22] = 30113; _outgoingV[22] = 323; _outgoingP[22] = -321;
_coupling[22] = 0.894/GeV; _maxweight[22] = 5.;
_incoming[23] = 30113; _outgoingV[23] = 313; _outgoingP[23] = -311;
_coupling[23] = 0.908/GeV; _maxweight[23] = 5.;
_incoming[24] = 30213; _outgoingV[24] = 323; _outgoingP[24] = -311;
_coupling[24] = 1.265/GeV; _maxweight[24] = 5.;
_incoming[25] = 30213; _outgoingV[25] = -313; _outgoingP[25] = 321;
_coupling[25] = 1.273/GeV; _maxweight[25] = 5.;
// omega'' rho pi
_incoming[26] = 30223; _outgoingV[26] = 213; _outgoingP[26] = -211;
_coupling[26] = 2.996/GeV; _maxweight[26] = 3.5;
_incoming[27] = 30223; _outgoingV[27] = 113; _outgoingP[27] = 111;
_coupling[27] = 2.996/GeV; _maxweight[27] = 3.5;
// omega' rho pi
_incoming[28] = 100223; _outgoingV[28] = 213; _outgoingP[28] = -211;
_coupling[28] = 4.507/GeV; _maxweight[28] = 4.;
_incoming[29] = 100223; _outgoingV[29] = 113; _outgoingP[29] = 111;
_coupling[29] = 4.507/GeV; _maxweight[29] = 4.;
// K*''->K* pi decays
_incoming[30] = 30313; _outgoingV[30] = 323; _outgoingP[30] = -211;
_coupling[30] = 3.36/GeV; _maxweight[30] = 4.;
_incoming[31] = 30313; _outgoingV[31] = 313; _outgoingP[31] = 111;
_coupling[31] = 2.38/GeV; _maxweight[31] = 4.;
_incoming[32] = 30323; _outgoingV[32] = 313; _outgoingP[32] = 211;
_coupling[32] = 3.36/GeV; _maxweight[32] = 4.;
_incoming[33] = 30323; _outgoingV[33] = 323; _outgoingP[33] = 111;
_coupling[33] = 2.38/GeV; _maxweight[33] = 4.;
// K*''->K rho decays
_incoming[34] = 30313; _outgoingP[34] = 321; _outgoingV[34] = -213;
_coupling[34] = 4.159/GeV; _maxweight[34] = 3.;
_incoming[35] = 30313; _outgoingP[35] = 311; _outgoingV[35] = 113;
_coupling[35] = 2.939/GeV; _maxweight[35] = 3.;
_incoming[36] = 30323; _outgoingP[36] = 311; _outgoingV[36] = 213;
_coupling[36] = 4.159/GeV; _maxweight[36] = 3.;
_incoming[37] = 30323; _outgoingP[37] = 321; _outgoingV[37] = 113;
_coupling[37] = 2.939/GeV; _maxweight[37] = 3.;
// K*' decays
_incoming[38] = 100313; _outgoingV[38] = 323; _outgoingP[38] = -211;
_coupling[38] = 9.469/GeV; _maxweight[38] = 6.;
_incoming[39] = 100313; _outgoingV[39] = 313; _outgoingP[39] = 111;
_coupling[39] = 6.781/GeV; _maxweight[39] = 6.;
_incoming[40] = 100323; _outgoingV[40] = 313; _outgoingP[40] = 211;
_coupling[40] = 9.469/GeV; _maxweight[40] = 6.;
_incoming[41] = 100323; _outgoingV[41] = 323; _outgoingP[41] = 111;
_coupling[41] = 6.781/GeV; _maxweight[41] = 6.;
// J/psi -> gamma eta_c decay
_incoming[42] = 443; _outgoingV[42] = 22; _outgoingP[42] = 441;
_coupling[42] = 0.149/GeV; _maxweight[42] = 30.;
// J/psi -> gamma eta' decay
_incoming[43] = 443; _outgoingV[43] = 22; _outgoingP[43] = 331;
_coupling[43] = 0.00250/GeV; _maxweight[43] = 1.7;
// J/psi -> rho pi decay
_incoming[44] = 443; _outgoingV[44] = 213; _outgoingP[44] = -211;
_coupling[44] = 0.00274/GeV; _maxweight[44] = 3.3;
_incoming[45] = 443; _outgoingV[45] = 113; _outgoingP[45] = 111;
_coupling[45] = 0.00274/GeV; _maxweight[45] = 3.3;
// J/psi -> K* K decay
_incoming[46] = 443; _outgoingV[46] = 323; _outgoingP[46] = -321;
_coupling[46] = 0.00180/GeV; _maxweight[46] = 6.;
_incoming[47] = 443; _outgoingV[47] = 313; _outgoingP[47] = -311;
_coupling[47] = 0.00180/GeV; _maxweight[47] = 6.;
// J/psi -> omega eta decay
_incoming[48] = 443; _outgoingV[48] = 223; _outgoingP[48] = 221;
_coupling[48] = 0.00154/GeV; _maxweight[48] = 6.;
// J/psi -> gamma eta decay
_incoming[49] = 443; _outgoingV[49] = 22; _outgoingP[49] = 221;
_coupling[49] = 0.00103/GeV; _maxweight[49] = 1.7;
// J/psi -> phi eta decay
_incoming[50] = 443; _outgoingV[50] = 333; _outgoingP[50] = 221;
_coupling[50] = 0.00110/GeV; _maxweight[50] = 5.5;
// J/psi -> phi eta' decay
_incoming[51] = 443; _outgoingV[51] = 333; _outgoingP[51] = 331;
_coupling[51] = 0.00085/GeV; _maxweight[51] = 5.5;
// J/psi -> omega pi0 decay
_incoming[52] = 443; _outgoingV[52] = 223; _outgoingP[52] = 111;
_coupling[52] = 0.00073/GeV; _maxweight[52] = 6.;
// J/psi -> rho0 eta decay
_incoming[53] = 443; _outgoingV[53] = 113; _outgoingP[53] = 221;
_coupling[53] = 0.00054/GeV; _maxweight[53] = 3.;
// J/psi -> rho0 eta' decay
_incoming[54] = 443; _outgoingV[54] = 113; _outgoingP[54] = 331;
_coupling[54] = 0.00045/GeV; _maxweight[54] = 3.;
// J/psi -> omega eta' decay
_incoming[55] = 443; _outgoingV[55] = 223; _outgoingP[55] = 331;
_coupling[55] = 0.00058/GeV; _maxweight[55] = 6.;
// J/psi -> gamma pi0 decay
_incoming[56] = 443; _outgoingV[56] = 22; _outgoingP[56] = 111;
_coupling[56] = 0.000177/GeV; _maxweight[56] = 1.7;
// psi(2s)-> j/psi eta decay
_incoming[57] = 100443; _outgoingV[57] = 443; _outgoingP[57] = 221;
_coupling[57] = 0.230/GeV; _maxweight[57] = 1.7;
// psi(2s)-> gamma eta_c decay
_incoming[58] = 100443; _outgoingV[58] = 22; _outgoingP[58] = 441;
_coupling[58] = 0.0114/GeV; _maxweight[58] = 3.4;
// psi(2s)-> gamma eta' decay
_incoming[59] = 100443; _outgoingV[59] = 22; _outgoingP[59] = 331;
_coupling[59] = 0.00062/GeV; _maxweight[59] = 1.7;
// psi(2s)-> J/psi pi0 decay
_incoming[60] = 100443; _outgoingV[60] = 443; _outgoingP[60] = 111;
_coupling[60] = 0.0106/GeV; _maxweight[60] = 1.7;
// psi(1d)-> gamma eta_c decay
_incoming[61] = 30443; _outgoingV[61] = 443; _outgoingP[61] = 221;
_coupling[61] = 0.135/GeV; _maxweight[61] = 1.7;
_incoming[62] = 30443; _outgoingV[62] = 333; _outgoingP[62] = 221;
_coupling[62] = 0.0076/GeV; _maxweight[62] = 5.5;
// psi(2s) K* K
_incoming[63] = 100443; _outgoingV[63] = 323; _outgoingP[63] = -321;
_coupling[63] = 0.00021/GeV; _maxweight[63] = 6.5;
_incoming[64] = 100443; _outgoingV[64] = 313; _outgoingP[64] = -311;
_coupling[64] = 0.00054/GeV; _maxweight[64] = 6.5;
// psi(2s) -> phi eta decay
_incoming[65] = 100443; _outgoingV[65] = 333; _outgoingP[65] = 221;
_coupling[65] = 0.00029/GeV; _maxweight[65] = 5.5;
// psi(2s) -> phi eta' decay
_incoming[66] = 100443; _outgoingV[66] = 333; _outgoingP[66] = 331;
_coupling[66] = 0.00033/GeV; _maxweight[66] = 5.5;
// psi(2s) -> rho pi decay
_incoming[67] = 100443; _outgoingV[67] = 213; _outgoingP[67] = -211;
_coupling[67] = 0.00017/GeV; _maxweight[67] = 3.5;
_incoming[68] = 100443; _outgoingV[68] = 113; _outgoingP[68] = 111;
_coupling[68] = 0.00017/GeV; _maxweight[68] = 3.5;
// psi(2s) -> omega eta' decay
_incoming[69] = 100443; _outgoingV[69] = 223; _outgoingP[69] = 331;
_coupling[69] = 0.00032/GeV; _maxweight[69] = 6.;
// psi(2s) -> rho0 eta decay
_incoming[70] = 100443; _outgoingV[70] = 113; _outgoingP[70] = 221;
_coupling[70] = 0.00025/GeV; _maxweight[70] = 3.5;
// psi(2s) -> omega pi0 decay
_incoming[71] = 100443; _outgoingV[71] = 223; _outgoingP[71] = 111;
_coupling[71] = 0.00022/GeV; _maxweight[71] = 6.;
// psi(2s) -> rho0 eta' decay
_incoming[72] = 100443; _outgoingV[72] = 113; _outgoingP[72] = 331;
_coupling[72] = 0.00026/GeV; _maxweight[72] = 3.5;
// initial size of the vectors for the database output
_initsize=_incoming.size();
}
int VectorMesonVectorPScalarDecayer::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(id ==_incoming[ix]) {
if((id1 ==_outgoingV[ix]&&id2 ==_outgoingP[ix])||
(id2 ==_outgoingV[ix]&&id1 ==_outgoingP[ix])) imode=ix;
}
if(idbar==_incoming[ix]) {
if((id1bar==_outgoingV[ix]&&id2bar==_outgoingP[ix])||
(id2bar==_outgoingV[ix]&&id1bar==_outgoingP[ix])) {
imode=ix;
cc=true;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
return imode;
}
void VectorMesonVectorPScalarDecayer::persistentOutput(PersistentOStream & os) const {
os << _incoming << _outgoingV << _outgoingP << _maxweight << ounit(_coupling,1/MeV);
}
void VectorMesonVectorPScalarDecayer::persistentInput(PersistentIStream & is, int) {
is >> _incoming >> _outgoingV >> _outgoingP >> _maxweight >> iunit(_coupling,1/MeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<VectorMesonVectorPScalarDecayer,DecayIntegrator>
describeHerwigVectorMesonVectorPScalarDecayer("Herwig::VectorMesonVectorPScalarDecayer", "HwVMDecay.so");
void VectorMesonVectorPScalarDecayer::Init() {
static ClassDocumentation<VectorMesonVectorPScalarDecayer> documentation
("The VectorMesonVectorPScalarDecayer class is designed for the "
"decay of a vector meson to another vector meson, or the photon, and a "
"pseudoscalar meson.");
static ParVector<VectorMesonVectorPScalarDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&VectorMesonVectorPScalarDecayer::_incoming,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMesonVectorPScalarDecayer,int> interfaceOutcomingVector
("OutgoingVector",
"The PDG code for the outgoing spin-1 particle",
&VectorMesonVectorPScalarDecayer::_outgoingV,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMesonVectorPScalarDecayer,int> interfaceOutcomingPScalar
("OutgoingPScalar",
"The PDG code for the outgoing spin-0 particle",
&VectorMesonVectorPScalarDecayer::_outgoingP,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<VectorMesonVectorPScalarDecayer,InvEnergy> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&VectorMesonVectorPScalarDecayer::_coupling,
1/MeV, 0, ZERO, -10000000/MeV, 10000000/MeV, false, false, true);
static ParVector<VectorMesonVectorPScalarDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&VectorMesonVectorPScalarDecayer::_maxweight,
0, 0, 0, -10000000, 10000000, false, false, true);
}
+
void VectorMesonVectorPScalarDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(_vectors[0],const_ptr_cast<tPPtr>(&part),
incoming,true,false);
VectorWaveFunction::constructSpinInfo(_vectors[1],decay[0],
outgoing,true,decay[0]->id()==ParticleID::gamma);
ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
}
+
double VectorMesonVectorPScalarDecayer::me2(const int,const Particle & part,
const tPDVector & ,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin1,PDT::Spin1,PDT::Spin0)));
// is the vector massless
bool photon(_outgoingV[imode()]==ParticleID::gamma);
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(_vectors[0],_rho,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
_vectors[1].resize(3);
for(unsigned int ix=0;ix<3;++ix) {
if(photon && ix==1) continue;
_vectors[1][ix] = HelicityFunctions::polarizationVector(-momenta[0],ix,Helicity::outgoing);
}
// compute the matrix element
LorentzPolarizationVector vtemp;
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1&&photon) {
for(unsigned int iy=0;iy<3;++iy) (*ME())(iy,ix,0)=0.;
}
else {
vtemp=_coupling[imode()]/part.mass()*
epsilon(part.momentum(),_vectors[1][ix],momenta[0]);
for(unsigned int iy=0;iy<3;++iy) (*ME())(iy,ix,0)=_vectors[0][iy].dot(vtemp);
}
}
double me = ME()->contract(_rho).real();
// test of the matrix element
// Energy pcm=Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(),
// momenta[1].mass());
// double test = sqr(_coupling[imode()]*pcm)*2./3.;
// cerr << "testing matrix element for " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << me << " " << test << " " << (me-test)/(me+test) << "\n";
// return the answer
return me;
}
bool VectorMesonVectorPScalarDecayer::twoBodyMEcode(const DecayMode & dm,int & mecode,
double & coupling) const {
int imode(-1);
int id(dm.parent()->id());
int idbar = dm.parent()->CC() ? dm.parent()->CC()->id() : id;
ParticleMSet::const_iterator pit(dm.products().begin());
int id1((**pit).id());
int id1bar = (**pit).CC() ? (**pit).CC()->id() : id1;
++pit;
int id2((**pit).id());
int id2bar = (**pit).CC() ? (**pit).CC()->id() : id2;
unsigned int ix(0); bool order(false);
do {
if(id==_incoming[ix]) {
if(id1 ==_outgoingV[ix]&&id2 ==_outgoingP[ix]) {
imode=ix;
order=true;
}
if(id2 ==_outgoingV[ix]&&id1 ==_outgoingP[ix]) {
imode=ix;
order=false;
}
}
if(idbar==_incoming[ix]&&imode<0) {
if(id1bar==_outgoingV[ix]&&id2bar==_outgoingP[ix]) {
imode=ix;
order=true;
}
if(id2bar==_outgoingV[ix]&&id1bar==_outgoingP[ix]) {
imode=ix;
order=false;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
coupling = _coupling[imode]*dm.parent()->mass();
mecode = 1;
return order;
}
// output the setup info for the particle database
void VectorMesonVectorPScalarDecayer::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() << ":OutgoingVector " << ix << " "
<< _outgoingV[ix] << "\n";
output << "newdef " << name() << ":OutgoingPScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*MeV << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "insert " << name() << ":OutgoingVector " << ix << " "
<< _outgoingV[ix] << "\n";
output << "insert " << name() << ":OutgoingPScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*MeV << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/VectorMeson/a1SimpleDecayer.cc b/Decay/VectorMeson/a1SimpleDecayer.cc
--- a/Decay/VectorMeson/a1SimpleDecayer.cc
+++ b/Decay/VectorMeson/a1SimpleDecayer.cc
@@ -1,453 +1,453 @@
// -*- C++ -*-
//
// a1SimpleDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the a1SimpleDecayer class.
//
#include "a1SimpleDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/WidthCalculatorBase.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
a1SimpleDecayer::a1SimpleDecayer()
// rho masses, widths and weights
: _rhomass ({0.773*GeV,1.370*GeV,1.750*GeV}),
_rhowidth({0.145*GeV,0.510*GeV,0.120*GeV}),
_rhowgts({1.0,-0.145,0.}),_localparameters(true),
_coupling(47.95/GeV),
// integration weights
_onemax(5.4474), _twomax(5.47784), _threemax(5.40185),
_onewgts ({0.235562,0.231098,0.131071,0.131135,0.135841,0.135294}),
_twowgts ({0.236208,0.229481,0.131169,0.133604,0.132685,0.136854}),
_threewgts({0.234259,0.233634,0.135922,0.129231,0.133949,0.133005}),
_mpi(ZERO) {
generateIntermediates(true);
}
void a1SimpleDecayer::doinit() {
DecayIntegrator::doinit();
// pointers to the particles we need as external particles
tPDPtr a1p = getParticleData(ParticleID::a_1plus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
tPDPtr pip = getParticleData(ParticleID::piplus);
tPDPtr pim = getParticleData(ParticleID::piminus);
tPDPtr pi0 = getParticleData(ParticleID::pi0);
// the different rho resonances
tPDPtr rhop[3] = {getParticleData(213),getParticleData(100213),
getParticleData(30213)};
tPDPtr rho0[3] = {getParticleData(113),getParticleData(100113),
getParticleData(30113)};
tPDPtr rhom[3] = {getParticleData(-213),getParticleData(-100213),
getParticleData(-30213)};
// decay mode a_1+ -> pi+ pi0 pi0
tPDPtr in = a1p;
tPDVector out = {pi0,pi0,pip};
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_onemax));
unsigned int nrho(0);
for(unsigned int ix=0;ix<3;++ix) if(rhop[ix]) ++nrho;
if(_onewgts.size()!=2*nrho) _onewgts=vector<double>(2*nrho,0.5/nrho);
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho+ channel
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rhop[ix],0,2,1,1,1,3));
c1.weight(_onewgts[2*ix]);
mode->addChannel(c1);
// second rho+ channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rhop[ix],0,1,1,2,1,3));
c2.weight(_onewgts[2*ix+1]);
mode->addChannel(c2);
}
addMode(mode);
// decay mode a_10 -> pi+ pi- pi0
in = a10;
out = {pip,pim,pi0};
mode = new_ptr(PhaseSpaceMode(in,out,_twomax));
if(_twowgts.size()!=2*nrho) _twowgts=vector<double>(2*nrho,0.5/nrho);
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho channel
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rhop[ix],0,2,1,1,1,3));
c1.weight(_twowgts[2*ix]);
mode->addChannel(c1);
// second channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rhom[ix],0,1,1,2,1,3));
c2.weight(_twowgts[2*ix+1]);
mode->addChannel(c2);
}
addMode(mode);
// decay mode a_1+ -> pi+ pi+ pi-
in = a1p;
out = {pip,pip,pim};
mode = new_ptr(PhaseSpaceMode(in,out,_threemax));
nrho = 0;
for(unsigned int ix=0;ix<3;++ix) if(rho0[ix]) ++nrho;
if(_threewgts.size()!=2*nrho) _threewgts=vector<double>(2*nrho,0.5/nrho);
for(unsigned int ix=0;ix<3;++ix) {
if(!rho0[ix]) continue;
// the neutral rho channels
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rho0[ix],0,2,1,1,1,3));
c1.weight(_threewgts[2*ix]);
mode->addChannel(c1);
// interchanged channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rho0[ix],0,1,1,2,1,3));
c2.weight(_threewgts[2*ix+1]);
mode->addChannel(c2);
}
addMode(mode);
// if using local parameters set the values in the phase space channels
if(_localparameters) {
for(unsigned int iy=0;iy<_rhomass.size();++iy) {
resetIntermediate(rho0[iy],_rhomass[iy],_rhowidth[iy]);
resetIntermediate(rhop[iy],_rhomass[iy],_rhowidth[iy]);
resetIntermediate(rhom[iy],_rhomass[iy],_rhowidth[iy]);
}
// make sure the rho array has enough masses
if(_rhomass.size()<3) {
for(unsigned int ix=_rhomass.size();ix<3;++ix) {
_rhomass.push_back(rhop[ix]->mass());
_rhowidth.push_back(rhop[ix]->width());
}
}
}
// set the local variables if needed
else {
// masses and widths for the particles
_rhomass.resize(3);_rhowidth.resize(3);
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
_rhomass[ix]=rhop[ix]->mass();
_rhowidth[ix]=rhop[ix]->width();
}
}
_mpi = pip->mass();
}
void a1SimpleDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
// get the weights for the different channels
for(unsigned int ix=0;ix<_onewgts.size();++ix)
_onewgts[ix]=mode(0)->channels()[ix].weight();
for(unsigned int ix=0;ix<_twowgts.size();++ix)
_twowgts[ix]=mode(1)->channels()[ix].weight();
for(unsigned int ix=0;ix<_threewgts.size();++ix)
_threewgts[ix]=mode(2)->channels()[ix].weight();
// get the maximum weight
_onemax = mode(0)->maxWeight();
_twomax = mode(1)->maxWeight();
_threemax = mode(2)->maxWeight();
}
}
void a1SimpleDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_rhomass,GeV) << ounit(_rhowidth,GeV) << _rhowgts
<< _localparameters << ounit(_coupling,1./GeV) << _onemax
<< _twomax << _threemax << _onewgts << _twowgts << _threewgts
<< ounit(_mpi,GeV);
}
void a1SimpleDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_rhomass,GeV) >> iunit(_rhowidth,GeV) >> _rhowgts
>> _localparameters >> iunit(_coupling,1./GeV) >> _onemax
>> _twomax >> _threemax >> _onewgts >> _twowgts >> _threewgts
>> iunit(_mpi,GeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<a1SimpleDecayer,DecayIntegrator>
describeHerwiga1SimpleDecayer("Herwig::a1SimpleDecayer", "HwVMDecay.so");
void a1SimpleDecayer::Init() {
static ClassDocumentation<a1SimpleDecayer> documentation
("The a1SimpleDecayer class implements a simple model for the decay of"
" the a_1 to three pions based on the approach of Kuhn and Santanmaria,"
" Z.Phys. C48, 445 (1990)",
"The decays of the $a_1$ were modelled using the approach of "
"\\cite{Kuhn:1990ad}.\n",
"\\bibitem{Kuhn:1990ad} J.~H.~Kuhn and A.~Santamaria,\n"
"Z.\\ Phys.\\ C {\\bf 48} (1990) 445.\n"
"%%CITATION = ZEPYA,C48,445;%%\n");
static Switch<a1SimpleDecayer,bool> interfaceLocalParameters
("LocalParameters",
"Use local values of the intermediate resonances masses and widths",
&a1SimpleDecayer::_localparameters, true, false, false);
static SwitchOption interfaceLocalParametersLocal
(interfaceLocalParameters,
"Local",
"Use the local values",
true);
static SwitchOption interfaceLocalParametersDefault
(interfaceLocalParameters,
"ParticleData",
"Use the values from the particleData objects",
false);
static Parameter<a1SimpleDecayer,InvEnergy> interfaceCoupling
("Coupling",
"The overall coupling for the decay",
&a1SimpleDecayer::_coupling, 1./GeV, 47.95/GeV, ZERO, 100./GeV,
false, false, Interface::limited);
static ParVector<a1SimpleDecayer,Energy> interfacerhomass
("RhoMasses",
"The masses of the different rho resonaces",
&a1SimpleDecayer::_rhomass, MeV, 3, 775.*MeV, ZERO, 10000*MeV,
false, false, Interface::limited);
static ParVector<a1SimpleDecayer,Energy> interfaceRhoWidths
("RhoWidths",
"The widths of the different rho resonances",
&a1SimpleDecayer::_rhowidth, MeV, 3, 141*MeV, 0.0*MeV, 10000.0*MeV,
false, false, Interface::limited);
static ParVector<a1SimpleDecayer,double> interfaceRhoWeights
("RhoWeights",
"Weight for the different rho resonances",
&a1SimpleDecayer::_rhowgts, 3, 0.0, -10.0, 10.0,
false, false, Interface::limited);
static Parameter<a1SimpleDecayer,double> interfaceOneMax
("OneMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi0pi0",
&a1SimpleDecayer::_onemax, 5.57613, 0.0, 10000.0,
false, false, true);
static Parameter<a1SimpleDecayer,double> interfaceTwoMax
("TwoMax",
"The maximum weight for the integration fo the channel a_1^0->pi+pi-pi0",
&a1SimpleDecayer::_twomax, 5.61218, 0.0, 10000.0,
false, false, true);
static Parameter<a1SimpleDecayer,double> interfaceThreeMax
("ThreeMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi+pi-",
&a1SimpleDecayer::_threemax, 5.5384, 0.0, 10000.0,
false, false, true);
static ParVector<a1SimpleDecayer,double> interfaceonewgts
("OneChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi0pi0",
&a1SimpleDecayer::_onewgts,
0, 0, 0, 0., 1., false, false, true);
static ParVector<a1SimpleDecayer,double> interfacetwowgts
("TwoChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^0->pi+pi-pi0",
&a1SimpleDecayer::_twowgts,
0, 0, 0, 0., 1., false, false, true);
static ParVector<a1SimpleDecayer,double> interfacethreewgts
("ThreeChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi+pi-",
&a1SimpleDecayer::_threewgts,
0, 0, 0, 0., 1., false, false, true);
}
int a1SimpleDecayer::modeNumber(bool & cc,tcPDPtr parent,
- const tPDVector & children) const {
+ const tPDVector & children) const {
if(children.size()!=3) return -1;
int id(parent->id());
// check the pions
int npi0(0),npiplus(0),npiminus(0);
for(auto child : children) {
int idtemp= child->id();
if(idtemp==ParticleID::piplus) ++npiplus;
else if(idtemp==ParticleID::piminus) ++npiminus;
else if(idtemp==ParticleID::pi0) ++npi0;
}
int imode(-1);
// a_1+ decay modes
if(id==ParticleID::a_1plus) {
cc=false;
if(npiplus==1&&npi0==2) imode=0;
else if(npiplus==2&&npiminus==1) imode=2;
}
// a_1- modes
else if(id==ParticleID::a_1minus) {
cc=true;
if(npiminus==1&&npi0==2) imode=0;
else if(npiminus==2&&npiplus==1) imode=2;
}
// a_0 modes
else if(id==ParticleID::a_10) {
cc=false;
if(npiminus==1&&npiplus==1&&npi0==1) imode=1;
}
return imode;
}
void a1SimpleDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(_vectors,const_ptr_cast<tPPtr>(&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 a1SimpleDecayer::me2(const int ichan, const Particle & part,
const tPDVector & ,
const vector<Lorentz5Momentum> & 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<tPPtr>(&part),
incoming,false);
}
Lorentz5Vector<complex<Energy> > current;
Energy2 s1 = (momenta[1]+momenta[2]).m2();
Energy2 s2 = (momenta[0]+momenta[2]).m2();
if(ichan<0) {
current = rhoFormFactor(s2,-1)*(momenta[0]-momenta[2])
+rhoFormFactor(s1,-1)*(momenta[1]-momenta[2]);
}
else if(ichan<3) {
current =
rhoFormFactor(s2,ichan)*(momenta[0]-momenta[2]);
}
else if(ichan<6) {
current =
rhoFormFactor(s1,-1)*(momenta[1]-momenta[2]);
}
// compute the matrix element
for(unsigned int ix=0;ix<3;++ix)
(*ME())(ix,0,0,0)=_coupling*current.dot(_vectors[ix]);
// matrix element and identical particle factor
double output=ME()->contract(_rho).real();
if(imode()!=1) output*=0.5;
// test the output
// Energy2 s3 = (momenta[0]+momenta[1]).m2();
// double test = threeBodyMatrixElement(imode(),sqr(part.mass()),
// s3,s2,s1,momenta[0].mass(),
// momenta[1].mass(),
// momenta[2].mass());
// if(ichan<0) cerr << "testing matrix element " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << outgoing[2]->PDGName() << output << " " << test << " "
// << (output-test)/(output+test) << "\n";
// return the answer
return output;
}
double a1SimpleDecayer::
threeBodyMatrixElement(const int iopt,const Energy2 q2, const Energy2 s3,
const Energy2 s2,const Energy2 s1, const Energy m1,
const Energy m2 ,const Energy m3) const {
Energy2 v12 = (s2-2.*sqr(m1)-2.*sqr(m3))+0.25*sqr(s1-s3-sqr(m1)+sqr(m3))/q2;
Energy2 v22 = (s1-2.*sqr(m2)-2.*sqr(m3))+0.25*sqr(s2-s3-sqr(m2)+sqr(m3))/q2;
Energy2 v1v2 = (0.5*q2-s3-0.5*(3*sqr(m3)-sqr(m1)-sqr(m2)))
+0.25*(s1-s3-sqr(m1)+sqr(m3))*(s2-s3-sqr(m2)+sqr(m3))/q2;
Complex rho1=rhoFormFactor(s2,-1);
Complex rho2=rhoFormFactor(s1,-1);
double me = sqr(_coupling)*real(v12*rho1*conj(rho1)+v22*rho2*conj(rho2)
+2.*v1v2*rho1*conj(rho2))/3.;
if(iopt!=1) me *= 0.5;
return me;
}
WidthCalculatorBasePtr
a1SimpleDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
int ncharged=0;
for( ; pit!=pend;++pit) {
if(abs((**pit).id())==ParticleID::piplus) ++ncharged;
}
--ncharged;
// integrator to perform the integral
vector<double> inweights;inweights.push_back(0.5);inweights.push_back(0.5);
vector<int> intype;intype.push_back(2);intype.push_back(3);
vector<Energy> inmass(2,_rhomass[0]),inwidth(2,_rhowidth[0]);
vector<double> inpow(2,0.0);
Energy mpi0=getParticleData(ParticleID::pi0)->mass();
Energy mpic=getParticleData(ParticleID::piplus)->mass();
Energy m[3];
m[0] = ncharged<2 ? mpi0 : mpic;
m[1] = m[0];
m[2] = (ncharged==0||ncharged==2) ? mpi0 : mpic;
return new_ptr(ThreeBodyAllOnCalculator<a1SimpleDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,ncharged,m[0],m[1],m[2]));
}
void a1SimpleDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
// parameters for the DecayIntegrator base class
DecayIntegrator::dataBaseOutput(output,false);
output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
output << "newdef " << name() << ":Coupling " << _coupling*GeV << "\n";
output << "newdef " << name() << ":OneMax " << _onemax << "\n";
output << "newdef " << name() << ":TwoMax " << _twomax << "\n";
output << "newdef " << name() << ":ThreeMax " << _threemax << "\n";
for(unsigned int ix=0;ix<_rhomass.size();++ix) {
output << "newdef " << name() << ":RhoMasses " << ix << " "
<< _rhomass[ix]/MeV << "\n";
}
for(unsigned int ix=0;ix<_rhowidth.size();++ix) {
output << "newdef " << name() << ":RhoWidths " << ix << " "
<< _rhowidth[ix]/MeV << "\n";
}
for(unsigned int ix=0;ix<_rhowgts.size();++ix) {
output << "newdef " << name() << ":RhoWeights " << ix << " "
<< _rhowgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_onewgts.size();++ix) {
output << "newdef " << name() << ":OneChargedWeights "
<< ix << " " << _onewgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_twowgts.size();++ix) {
output << "newdef " << name() << ":TwoChargedWeights "
<< ix << " " << _twowgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_threewgts.size();++ix) {
output << "newdef " << name() << ":ThreeChargedWeights "
<< ix << " " << _threewgts[ix] << "\n";
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
// functions to return the Breit-Wigners
Complex a1SimpleDecayer::rhoFormFactor(Energy2 q2,int ires) const {
Complex output(0.),norm(0.);
for(unsigned int ix=0;ix<3;++ix) norm += _rhowgts[ix];
if(ires<0) {
for(unsigned int ix=0;ix<3;++ix) {
output+=_rhowgts[ix]*rhoBreitWigner(q2,ix);
}
}
else {
assert(ires<3);
output=_rhowgts[ires]*rhoBreitWigner(q2,ires);
}
return output/norm;
}
diff --git a/Decay/VectorMeson/a1ThreePionDecayer.cc b/Decay/VectorMeson/a1ThreePionDecayer.cc
--- a/Decay/VectorMeson/a1ThreePionDecayer.cc
+++ b/Decay/VectorMeson/a1ThreePionDecayer.cc
@@ -1,757 +1,735 @@
// -*- C++ -*-
//
// a1ThreePionDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the a1ThreePionDecayer class.
//
#include "a1ThreePionDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.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;
using namespace ThePEG::Helicity;
a1ThreePionDecayer::a1ThreePionDecayer()
: _rhomass(1,0.7761*GeV), _rhowidth(1,0.1445*GeV), _sigmamass(0.8*GeV),
_sigmawidth(0.8*GeV), _psigma(ZERO), _mpi(ZERO), _mpi2(ZERO),
_lambda2(1.2*GeV2), _a1mass2(1.23*1.23*GeV2),
_zsigma(0.), _zmag(1.3998721), _zphase(0.43585036),
_rhomag(1,1.), _rhophase(1,0.), _coupling(90.44),
- _localparameters(true), _zerowgts(3), _onewgts(7), _twowgts(7),
- _threewgts(8) ,_zeromax(19.144), _onemax(7.83592),
+ _localparameters(true),
+ _zerowgts ({0.339108,0.335601,0.325291}),
+ _onewgts ({0.19616 ,0.191408,0.12137 ,0.115498,0.12729 ,0.127183,0.12109 }),
+ _twowgts ({0.188163,0.192479,0.121658,0.12135 ,0.127298,0.124835,0.124217}),
+ _threewgts({0.153071,0.165741,0.107509,0.10275 ,0.109738,0.11254 ,0.125344,0.123307}) ,
+ _zeromax(19.144), _onemax(7.83592),
_twomax(6.64804), _threemax(6.66296) {
- // weights for the different channels
- _threewgts[0] = 0.153071;
- _threewgts[1] = 0.165741;
- _threewgts[2] = 0.107509;
- _threewgts[3] = 0.10275 ;
- _threewgts[4] = 0.109738;
- _threewgts[5] = 0.11254 ;
- _threewgts[6] = 0.125344;
- _threewgts[7] = 0.123307;
- _onewgts[0] = 0.19616 ;
- _onewgts[1] = 0.191408;
- _onewgts[2] = 0.12137 ;
- _onewgts[3] = 0.115498;
- _onewgts[4] = 0.12729 ;
- _onewgts[5] = 0.127183;
- _onewgts[6] = 0.12109 ;
- _zerowgts[0] = 0.339108;
- _zerowgts[1] = 0.335601;
- _zerowgts[2] = 0.325291;
- _twowgts[0] = 0.188163;
- _twowgts[1] = 0.192479;
- _twowgts[2] = 0.121658;
- _twowgts[3] = 0.12135 ;
- _twowgts[4] = 0.127298;
- _twowgts[5] = 0.124835;
- _twowgts[6] = 0.124217;
// generation of intermediates
generateIntermediates(true);
}
void a1ThreePionDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
// get the weights for the different channels
for(unsigned int ix=0;ix<_zerowgts.size();++ix)
_zerowgts[ix]=mode(0)->channels()[ix].weight();
for(unsigned int ix=0;ix<_onewgts.size();++ix)
_onewgts[ix]=mode(1)->channels()[ix].weight();
for(unsigned int ix=0;ix<_twowgts.size();++ix)
_twowgts[ix]=mode(2)->channels()[ix].weight();
for(unsigned int ix=0;ix<_threewgts.size();++ix)
_threewgts[ix]=mode(3)->channels()[ix].weight();
// get the maximum weight
_zeromax = mode(0)->maxWeight();
_onemax = mode(1)->maxWeight();
_twomax = mode(2)->maxWeight();
_threemax = mode(3)->maxWeight();
}
}
void a1ThreePionDecayer::doinit() {
DecayIntegrator::doinit();
// particles we need for the external state
tPDPtr a1p = getParticleData(ParticleID::a_1plus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
tPDPtr pip = getParticleData(ParticleID::piplus);
tPDPtr pim = getParticleData(ParticleID::piminus);
tPDPtr pi0 = getParticleData(ParticleID::pi0);
// possible intermediate particles
// the different rho resonances
tPDPtr rhop[3] = {getParticleData(213),getParticleData(100213),
getParticleData(30213)};
tPDPtr rho0[3] = {getParticleData(113),getParticleData(100113),
getParticleData(30113)};
tPDPtr rhom[3] = {getParticleData(-213),getParticleData(-100213),
getParticleData(-30213)};
// the sigma
tPDPtr sigma = getParticleData(9000221);
// set up the phase space integration
// decay mode a_0 -> pi0 pi0 pi0
tPDPtr in = a10;
tPDVector out={pi0,pi0,pi0};
if(sigma) {
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_zeromax));
if(_zerowgts.size()!=3) _zerowgts=vector<double>(3,1./3.);
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,sigma,0,1,1,2,1,3));
c1.weight(_zerowgts[0]);
mode->addChannel(c1);
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,sigma,0,2,1,1,1,3));
c2.weight(_zerowgts[1]);
mode->addChannel(c2);
PhaseSpaceChannel c3((PhaseSpaceChannel(mode),0,sigma,0,3,1,1,1,2));
c3.weight(_zerowgts[2]);
mode->addChannel(c3);
addMode(mode);
}
// decay mode a_1+ -> pi+ pi0 pi0
in = a1p;
out = {pi0,pi0,pip};
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_onemax));
unsigned int nrho(0);
for(unsigned int ix=0;ix<3;++ix) if(rhop[ix]) ++nrho;
if(_onewgts.size()!=2*nrho+1) _onewgts=vector<double>(2*nrho+1,1./(2.*nrho+1.));
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho+ channel
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rhop[ix],0,1,1,2,1,3));
c1.weight(_onewgts[2*ix]);
mode->addChannel(c1);
// second rho+ channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rhop[ix],0,2,1,1,1,3));
c2.weight(_onewgts[2*ix+1]);
mode->addChannel(c2);
}
// the sigma channel
if(sigma) {
PhaseSpaceChannel c3((PhaseSpaceChannel(mode),0,sigma,0,3,1,1,1,2));
c3.weight(_onewgts.back());
mode->addChannel(c3);
}
addMode(mode);
// decay mode a_1 -> pi+ pi- pi0
in = a10;
out = {pip,pim,pi0};
mode = new_ptr(PhaseSpaceMode(in,out,_twomax));
if(_twowgts.size()!=2*nrho+1) _twowgts=vector<double>(2*nrho+1,1./(2.*nrho+1.));
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho channel
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rhop[ix],0,1,1,2,1,3));
c1.weight(_twowgts[2*ix]);
mode->addChannel(c1);
// second channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rhom[ix],0,2,1,1,1,3));
c2.weight(_twowgts[2*ix+1]);
mode->addChannel(c2);
}
// sigma channel
if(sigma) {
PhaseSpaceChannel c3((PhaseSpaceChannel(mode),0,sigma,0,3,1,1,1,2));
c3.weight(_twowgts.back());
mode->addChannel(c3);
}
addMode(mode);
// decay mode a_1+ -> pi+ pi+ pi-
in = a1p;
out = {pip,pip,pim};
mode = new_ptr(PhaseSpaceMode(in,out,_threemax));
nrho = 0;
for(unsigned int ix=0;ix<3;++ix) if(rho0[ix]) ++nrho;
if(_threewgts.size()!=2*nrho+2) _threewgts=vector<double>(2*nrho+2,1./(2.*nrho+2.));
for(unsigned int ix=0;ix<3;++ix) {
// the neutral rho channels
if(!rho0[ix]) continue;
// the neutral rho channels
PhaseSpaceChannel c1((PhaseSpaceChannel(mode),0,rho0[ix],0,1,1,2,1,3));
c1.weight(_threewgts[2*ix]);
mode->addChannel(c1);
// interchanged channel
PhaseSpaceChannel c2((PhaseSpaceChannel(mode),0,rho0[ix],0,2,1,1,1,3));
c2.weight(_threewgts[2*ix+1]);
mode->addChannel(c2);
}
// the sigma channels
if(sigma) {
PhaseSpaceChannel c3((PhaseSpaceChannel(mode),0,sigma,0,1,1,2,1,3));
c3.weight(_threewgts[6]);
mode->addChannel(c3);
PhaseSpaceChannel c4((PhaseSpaceChannel(mode),0,sigma,0,2,1,1,1,3));
c4.weight(_threewgts[7]);
mode->addChannel(c4);
}
addMode(mode);
// set up the parameters
_mpi=getParticleData(ParticleID::piplus)->mass();
_mpi2=sqr(_mpi);
if(_localparameters) {
if(_rhomass.size()<_rhocoupling.size()) {
unsigned int itemp=_rhomass.size();
_rhomass.resize(_rhocoupling.size());
_rhowidth.resize(_rhocoupling.size());
for(unsigned int ix=itemp;ix<_rhocoupling.size();++ix) {
_rhomass[ix]=rhop[ix]->mass();
_rhowidth[ix]=rhop[ix]->width();
}
// reset the intermediates in the phase space integration if needed
resetIntermediate(sigma,_sigmamass,_sigmawidth);
for(unsigned int iy=0;iy<_rhocoupling.size();++iy) {
resetIntermediate(rho0[iy],_rhomass[iy],_rhowidth[iy]);
resetIntermediate(rhop[iy],_rhomass[iy],_rhowidth[iy]);
resetIntermediate(rhom[iy],_rhomass[iy],_rhowidth[iy]);
}
}
}
else {
_a1mass2=sqr(getParticleData(ParticleID::a_1plus)->mass());
if(sigma) {
_sigmamass=sigma->mass();
_sigmawidth=sigma->width();
}
_rhomass.resize(_rhocoupling.size());
_rhowidth.resize(_rhocoupling.size());
for(unsigned int ix=0;ix<_rhocoupling.size();++ix) {
_rhomass[ix]=rhop[ix]->mass();
_rhowidth[ix]=rhop[ix]->width();
}
}
// parameters for the resonances
// for the sigma
_psigma=Kinematics::pstarTwoBodyDecay(_sigmamass,_mpi,_mpi);
// for the rho
_prho.resize(_rhomass.size());_hm2.resize(_rhomass.size());
_dhdq2m2.resize(_rhomass.size());_rhoD.resize(_rhomass.size());
for(unsigned int ix=0;ix<_rhomass.size();++ix) {
_prho[ix] = Kinematics::pstarTwoBodyDecay(_rhomass[ix],_mpi,_mpi);
_hm2[ix] = hFunction(_rhomass[ix]);
_dhdq2m2[ix] = dhdq2Parameter(ix);
_rhoD[ix] = DParameter(ix);
}
// convert the magnitude and phase of z into a phase
_zsigma = _zmag*Complex(cos(_zphase),sin(_zphase));
// convert rho couplings
for(unsigned int ix=0;ix<_rhomag.size();++ix) {
_rhocoupling.push_back(_rhomag[ix]*Complex(cos(_rhophase[ix]),sin(_rhophase[ix])));
}
}
int a1ThreePionDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
if(children.size()!=3) return -1;
int id(parent->id());
// check the pions
tPDVector::const_iterator pit = children.begin();
tPDVector::const_iterator pend = children.end();
int idtemp,npi0(0),npiplus(0),npiminus(0);
for( ; pit!=pend;++pit) {
idtemp=(**pit).id();
if(idtemp==ParticleID::piplus) ++npiplus;
else if(idtemp==ParticleID::piminus) ++npiminus;
else if(idtemp==ParticleID::pi0) ++npi0;
}
int imode(-1);
// a_1+ decay modes
if(id==ParticleID::a_1plus) {
cc=false;
if(npiplus==1&&npi0==2) imode=1;
else if(npiplus==2&&npiminus==1) imode=3;
}
// a_1- modes
else if(id==ParticleID::a_1minus) {
cc=true;
if(npiminus==1&&npi0==2) imode=1;
else if(npiminus==2&&npiplus==1) imode=3;
}
// a_0 modes
else if(id==ParticleID::a_10) {
cc=false;
if(npiminus==1&&npiplus==1&&npi0==1) imode=2;
else if(npi0==3) imode=0;
}
return imode;
}
void a1ThreePionDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_rhomass,GeV) << ounit(_rhowidth,GeV) << ounit(_prho,GeV)
<< ounit(_hm2,GeV2) << ounit(_rhoD,GeV2) << _dhdq2m2 << ounit(_sigmamass,GeV)
<< ounit(_sigmawidth,GeV) << ounit(_psigma,GeV) << ounit(_mpi,GeV)
<< ounit(_mpi2,GeV2) << ounit(_lambda2,GeV2) << ounit(_a1mass2,GeV2) << _zsigma
<< _rhocoupling << _coupling << _localparameters << _zerowgts << _onewgts
<< _twowgts << _threewgts << _zeromax << _zmag << _zphase
<< _onemax << _twomax << _threemax << _coupling << _rhomag << _rhophase;
}
void a1ThreePionDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_rhomass,GeV) >> iunit(_rhowidth,GeV) >> iunit(_prho,GeV)
>> iunit(_hm2,GeV2) >> iunit(_rhoD,GeV2) >> _dhdq2m2 >> iunit(_sigmamass,GeV)
>> iunit(_sigmawidth,GeV) >> iunit(_psigma,GeV) >> iunit(_mpi,GeV)
>> iunit(_mpi2,GeV2) >> iunit(_lambda2,GeV2) >> iunit(_a1mass2,GeV2) >> _zsigma
>> _rhocoupling >> _coupling >> _localparameters >> _zerowgts >> _onewgts
>> _twowgts >> _threewgts >> _zeromax >> _zmag >> _zphase
>> _onemax >> _twomax >> _threemax >> _coupling >> _rhomag >> _rhophase;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<a1ThreePionDecayer,DecayIntegrator>
describeHerwiga1ThreePionDecayer("Herwig::a1ThreePionDecayer", "HwVMDecay.so");
void a1ThreePionDecayer::Init() {
static ClassDocumentation<a1ThreePionDecayer> documentation
("The a1ThreePionDecayer class is designed to decay the a_1 "
"resonance to three pions using a model based on that used in the modelling "
"of tau->4 pions.",
"The decay of the $a_1$ resonance to three pions uses a model based on"
"tau to four pions, \\cite{Bondar:2002mw}.",
"%\\cite{Bondar:2002mw}\n"
"\\bibitem{Bondar:2002mw}\n"
" A.~E.~Bondar, S.~I.~Eidelman, A.~I.~Milstein, T.~Pierzchala, N.~I.~Root, Z.~Was and M.~Worek,\n"
" ``Novosibirsk hadronic currents for tau --> 4pi channels of tau decay\n"
" %library TAUOLA,''\n"
" Comput.\\ Phys.\\ Commun.\\ {\\bf 146}, 139 (2002)\n"
" [arXiv:hep-ph/0201149].\n"
" %%CITATION = CPHCB,146,139;%%\n"
);
static Switch<a1ThreePionDecayer,bool> interfaceLocalParameters
("LocalParameters",
"Use local values of the intermediate resonances masses and widths",
&a1ThreePionDecayer::_localparameters, true, false, false);
static SwitchOption interfaceLocalParametersLocal
(interfaceLocalParameters,
"Local",
"Use the local values",
true);
static SwitchOption interfaceLocalParametersDefault
(interfaceLocalParameters,
"ParticleData",
"Use the values from the particleData objects",
false);
static Parameter<a1ThreePionDecayer,double> interfaceCoupling
("Coupling",
"The overall coupling for the decay",
&a1ThreePionDecayer::_coupling, 90.44, 0.0, 1000.0,
false, false, true);
static ParVector<a1ThreePionDecayer,Energy> interfacerhomass
("RhoMasses",
"The masses of the different rho resonnaces",
&a1ThreePionDecayer::_rhomass,
GeV, 0, ZERO, ZERO, 10000*GeV, false, false, true);
static ParVector<a1ThreePionDecayer,Energy> interfacerhowidth
("RhoWidths",
"The widths of the different rho resonnaces",
&a1ThreePionDecayer::_rhowidth,
GeV, 0, ZERO, ZERO, 10000*GeV, false, false, true);
static ParVector<a1ThreePionDecayer,double> interfaceRhoMagnitude
("RhoMagnitude",
"The magnitude of the rho couplings",
&a1ThreePionDecayer::_rhomag, -1, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static ParVector<a1ThreePionDecayer,double> interfaceRhoPhase
("RhoPhase",
"The phase of the rho coupling",
&a1ThreePionDecayer::_rhophase, -1, 0., 0.0, 2.*Constants::pi,
false, false, Interface::limited);
static Parameter<a1ThreePionDecayer,Energy2> interfaceLambda2
("Lambda2",
"The value of the mass scale squared to use in the form-factor",
&a1ThreePionDecayer::_lambda2, GeV2, 1.2*GeV2, 0.0001*GeV2, 10.0*GeV2,
false, false, true);
static Parameter<a1ThreePionDecayer,Energy2> interfacea1mass2
("a1mass2",
"The local value of the square of the a_1 mass",
&a1ThreePionDecayer::_a1mass2, GeV2, 1.5129*GeV2, 0.5*GeV2, 10.0*GeV2,
false, false, true);
static Parameter<a1ThreePionDecayer,Energy> interfaceSigmaMass
("SigmaMass",
"The local value of the sigma mass",
&a1ThreePionDecayer::_sigmamass, GeV, 0.8*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<a1ThreePionDecayer,Energy> interfaceSigmaWidth
("SigmaWidth",
"The local value of the sigma width",
&a1ThreePionDecayer::_sigmawidth, GeV, 0.8*GeV, ZERO, 10.0*GeV,
false, false, true);
static ParVector<a1ThreePionDecayer,double> interfacezerowgts
("AllNeutralWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^0->pi0pi0pi0",
&a1ThreePionDecayer::_zerowgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionDecayer,double> interfaceonewgts
("OneChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi0pi0",
&a1ThreePionDecayer::_onewgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionDecayer,double> interfacetwowgts
("TwoChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^0->pi+pi-pi0",
&a1ThreePionDecayer::_twowgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionDecayer,double> interfacethreewgts
("ThreeChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi+pi-",
&a1ThreePionDecayer::_threewgts,
0, 0, 0, -10000, 10000, false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceZeroMax
("ZeroMax",
"The maximum weight for the integration fo the channel a_1^0->pi0pi0pi0",
&a1ThreePionDecayer::_zeromax, 107.793, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceOneMax
("OneMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi0pi0",
&a1ThreePionDecayer::_onemax, 1088.96, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceTwoMax
("TwoMax",
"The maximum weight for the integration fo the channel a_1^0->pi+pi-pi0",
&a1ThreePionDecayer::_twomax, 1750.73, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceThreeMax
("ThreeMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi+pi-",
&a1ThreePionDecayer::_threemax, 739.334, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceSigmaMagnitude
("SigmaMagnitude",
"magnitude of the relative sigma coupling",
&a1ThreePionDecayer::_zmag, 1.3998721, 0.0, 10.0e20,
false, false, true);
static Parameter<a1ThreePionDecayer,double> interfaceSigmaPhase
("SigmaPhase",
"phase of the relative sigma coupling",
&a1ThreePionDecayer::_zphase, 0.43585036, 0.0, Constants::twopi,
false, false, true);
}
void a1ThreePionDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(_vectors,const_ptr_cast<tPPtr>(&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 a1ThreePionDecayer::me2(const int ichan, const Particle & part,
const tPDVector & ,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
useMe();
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin0,PDT::Spin0,PDT::Spin0)));
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(_vectors,_rho,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
// momentum of the incoming particle
Lorentz5Momentum Q=part.momentum();
// momenta of the intermediates
Energy2 s1=(momenta[1]+momenta[2]).m2();
Energy2 s2=(momenta[0]+momenta[2]).m2();
Energy2 s3=(momenta[0]+momenta[1]).m2();
Energy2 dot01=Q*momenta[0];
Energy2 dot02=Q*momenta[1];
Energy2 dot03=Q*momenta[2];
// vector for the output
LorentzVector<complex<Energy3> > output;
// a_10 -> pi0pi0pi0
if(imode()==0) {
//the breit-wigners
Complex sig1=sigmaBreitWigner(s1);
Complex sig2=sigmaBreitWigner(s2);
Complex sig3=sigmaBreitWigner(s3);
// compute the vector
LorentzPolarizationVectorE tmpoutput;
if(ichan<0) {
tmpoutput= sig1*(momenta[0])+sig2*(momenta[1])
+sig3*(momenta[2]);
}
else if(ichan==0) tmpoutput=sig1*(momenta[0]);
else if(ichan==1) tmpoutput=sig2*(momenta[1]);
else if(ichan==2) tmpoutput=sig3*(momenta[2]);
// the coupling z and identical particle factor
output = tmpoutput * _zsigma* 1./sqrt(6.) *Q.mass2();
}
// a_1+ -> pi0pi0pi+
else if(imode()==1) {
// scalar propagator
Complex sig1 = sigmaBreitWigner(s3);
// sigma terms
if(ichan<0||ichan==6)
output = _zsigma*Q.mass2()*sig1*momenta[2];
// the rho terms
for(int ix=0,N=_rhocoupling.size();ix<N;++ix) {
Complex rho1=_rhocoupling[ix]*rhoBreitWigner(s1,ix);
Complex rho2=_rhocoupling[ix]*rhoBreitWigner(s2,ix);
if(ichan<0||ichan==2*ix) {
output +=rho1*(dot03*(momenta[1])-
dot02*(momenta[2]));
}
if(ichan<0||ichan==2*ix+1){
output +=rho2*(dot03*(momenta[0])-
dot01*(momenta[2]));
}
}
// the identical particle factor
output *= 1./sqrt(2.);
}
// a_10->pi+pi-pi0
else if(imode()==2) {
// the sigma terms
Complex sig1=sigmaBreitWigner(s3);
if(ichan<0||ichan==6)
output = _zsigma*Q.mass2()*sig1*momenta[2];
// rho terms
for(int ix=0,N=_rhocoupling.size();ix<N;++ix) {
Complex rho1=_rhocoupling[ix]*rhoBreitWigner(s1,ix);
Complex rho2=_rhocoupling[ix]*rhoBreitWigner(s2,ix);
if(ichan<0||ichan==2*ix) {
output+=rho1*(dot03*(momenta[1])
-dot02*(momenta[2]));
}
if(ichan<0||ichan==2*ix+1) {
output+=rho2*(dot03*(momenta[0])
-dot01*(momenta[2]));
}
}
}
// a1+ -> pi+pi+pi-
else if(imode()==3) {
// the scalar propagators
Complex sig1=sigmaBreitWigner(s1);
Complex sig2=sigmaBreitWigner(s2);
// sigma terms
LorentzPolarizationVectorE tmpoutput;
if(ichan<0||ichan==6) tmpoutput+=sig1*(momenta[0]);
if(ichan<0||ichan==7) tmpoutput+=sig2*(momenta[1]);
output = tmpoutput * _zsigma * Q.mass2();
// rho terms
for(int ix=0,N=_rhocoupling.size();ix<N;++ix) {
Complex rho1 = _rhocoupling[ix]*rhoBreitWigner(s1,ix);
Complex rho2 = _rhocoupling[ix]*rhoBreitWigner(s2,ix);
if(ichan<0||ichan==2*ix) {
output-=rho1*( dot03*(momenta[1])-
dot02*(momenta[2]));
}
if(ichan<0||ichan==2*ix+1) {
output-=rho2*( dot03*(momenta[0])-
dot01*(momenta[2]));
}
}
// the identical particle factor
output *= 1./sqrt(2.);
}
// form-factor
LorentzPolarizationVector outputFinal
= output * a1FormFactor(Q.mass2())*_coupling/(Q.mass()*sqr(_rhomass[0]));
// compute the matrix element
for(unsigned int ix=0;ix<3;++ix)
(*ME())(ix,0,0,0)=outputFinal.dot(_vectors[ix]);
// return the answer
double out = ME()->contract(_rho).real();
// test of the answer
// double test = threeBodyMatrixElement(imode(),sqr(part.mass()),s3,s2,s1,
// momenta[0].mass(),momenta[1].mass(),
// momenta[2].mass());
// if(ichan<0) cerr << "testing matrix element " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << outgoing[2]->PDGName() << " " << out << " " << test << " "
// << (out-test)/(out+test) << "\n";
return out;
}
// matrix element for the running a_1 width
double a1ThreePionDecayer::
threeBodyMatrixElement(const int iopt,const Energy2 q2, const Energy2 s3,
const Energy2 s2,const Energy2 s1, const Energy m1,
const Energy m2 ,const Energy m3) const {
Energy6 meout(0.*pow<3,1>(GeV2));
Energy2 m12(sqr(m1)),m22(sqr(m2)),m32(sqr(m3));
Energy2 dot01(q2-s1+m12),dot02(q2-s2+m22),dot03(q2-s3+m32),
dot12(s3-m12-m22),dot13(s2-m12-m32),dot23(s1-m22-m32);
if(iopt==0) {
Complex sig1=sigmaBreitWigner(s1);
Complex sig2=sigmaBreitWigner(s2);
Complex sig3=sigmaBreitWigner(s3);
Energy2 metemp =
real(0.25*sig1*conj(sig1)*lambda(q2,s1,m12)/q2+
0.25*sig2*conj(sig2)*lambda(q2,s2,m22)/q2+
0.25*sig3*conj(sig3)*lambda(q2,s3,m32)/q2+
sig1*conj(sig2)*(-dot12+0.5*dot01*dot02/q2)+
sig1*conj(sig3)*(-dot13+0.5*dot01*dot03/q2)+
sig2*conj(sig3)*(-dot23+0.5*dot02*dot03/q2));
meout = metemp*real(_zsigma*conj(_zsigma))/6.*sqr(q2);
}
else if(iopt==1||iopt==2) {
// the sigma terms
Complex sig=sigmaBreitWigner(s3);
Complex rho1,rho2;
for(int ix=0,N=_rhocoupling.size();ix<N;++ix) {
rho1 += _rhocoupling[ix]*rhoBreitWigner(s1,ix);
rho2 += _rhocoupling[ix]*rhoBreitWigner(s2,ix);
}
meout =
0.25*lambda(q2,m32,s3)*q2*norm(_zsigma*sig)+
0.25*norm(rho1)*(dot23*dot02*dot03-m32*sqr(dot02)-m22*sqr(dot03))+
0.25*norm(rho2)*(dot13*dot01*dot03-m32*sqr(dot01)-m12*sqr(dot03))-
0.5*real(_zsigma*sig*conj(rho1))*q2*(dot03*dot23-2.*m32*dot02)-
0.5*real(_zsigma*sig*conj(rho2))*q2*(dot03*dot13-2.*m32*dot01)-
0.25*real(rho1*conj(rho2))*(sqr(dot03)*dot12-dot03*dot02*dot13
-dot03*dot01*dot23+2.*m32*dot02*dot01);
if(iopt==1) meout *= 0.5;
}
else if(iopt==3) {
Complex sig1=sigmaBreitWigner(s1);
Complex sig2=sigmaBreitWigner(s2);
Complex rho1,rho2;
for(int ix=0,N=_rhocoupling.size();ix<N;++ix) {
rho1 += _rhocoupling[ix]*rhoBreitWigner(s1,ix);
rho2 += _rhocoupling[ix]*rhoBreitWigner(s2,ix);
}
meout =
0.25*lambda(q2,m12,s1)*q2*norm(_zsigma*sig1)+
0.25*lambda(q2,m22,s2)*q2*norm(_zsigma*sig2)+
0.25*norm(rho1)*(dot23*dot02*dot03-m32*sqr(dot02)-m22*sqr(dot03))+
0.25*norm(rho2)*(dot13*dot01*dot03-m32*sqr(dot01)-m12*sqr(dot03))-
0.25*real(rho1*conj(rho2))*(sqr(dot03)*dot12-dot03*dot02*dot13
-dot03*dot01*dot23+2.*m32*dot02*dot01)-
real(_zsigma*sig1*conj(_zsigma*sig2))*q2*(q2*dot12-0.5*dot02*dot01)+
0.5*real(_zsigma*sig1*conj(rho1))*q2*(dot03*dot12-dot02*dot13)+
0.5*real(_zsigma*sig2*conj(rho1))*q2*(2.*dot03*m22-dot02*dot23)+
0.5*real(_zsigma*sig1*conj(rho2))*q2*(2.*dot03*m12-dot01*dot13)+
0.5*real(_zsigma*sig2*conj(rho2))*q2*(dot03*dot12-dot01*dot23);
meout *= 0.5;
}
return meout*a1FormFactor(q2)*sqr(_coupling/sqr(_rhomass[0]))/q2/3.;
}
WidthCalculatorBasePtr
a1ThreePionDecayer::threeBodyMEIntegrator(const DecayMode & dm) const {
ParticleMSet::const_iterator pit = dm.products().begin();
ParticleMSet::const_iterator pend = dm.products().end();
int ncharged=0;
for( ; pit!=pend;++pit) {
if(abs((**pit).id())==ParticleID::piplus) ++ncharged;
}
// integrator to perform the integral
vector<double> inweights(2,0.5);
vector<int> intype;intype.push_back(2);intype.push_back(3);
vector<Energy> inmass(2,_rhomass[0]),inwidth(2,_rhowidth[0]);
vector<double> inpow(2,0.0);
Energy m[3];
Energy mpi0=getParticleData(ParticleID::pi0)->mass();
Energy mpic=getParticleData(ParticleID::piplus)->mass();
m[0] = ncharged<2 ? mpi0 : mpic;
m[1] = m[0];
m[2] = (ncharged==0||ncharged==2) ? mpi0 : mpic;
return new_ptr(ThreeBodyAllOnCalculator<a1ThreePionDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,ncharged,m[0],m[1],m[2]));
}
// output the setup information for the particle database
void a1ThreePionDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
// parameters for the DecayIntegrator base class
DecayIntegrator::dataBaseOutput(output,false);
output << "newdef " << name() << ":LocalParameters " << _localparameters << "\n";
output << "newdef " << name() << ":Coupling " << _coupling << "\n";
output << "newdef " << name() << ":Lambda2 " << _lambda2/GeV2 << "\n";
output << "newdef " << name() << ":a1mass2 " << _a1mass2/GeV2 << "\n";
output << "newdef " << name() << ":SigmaMass " << _sigmamass/GeV << "\n";
output << "newdef " << name() << ":SigmaWidth " << _sigmawidth/GeV << "\n";
output << "newdef " << name() << ":SigmaMagnitude " << _zmag << "\n";
output << "newdef " << name() << ":SigmaPhase " << _zphase << "\n";
for(unsigned int ix=0;ix<_rhomag.size();++ix) {
if(ix<1) output << "newdef " << name() << ":RhoMagnitude " << ix << " "
<< _rhomag[ix] << "\n";
else output << "insert " << name() << ":RhoMagnitude " << ix << " "
<< _rhomag[ix] << "\n";
}
for(unsigned int ix=0;ix<_rhophase.size();++ix) {
if(ix<1) output << "newdef " << name() << ":RhoPhase " << ix << " "
<< _rhophase[ix] << "\n";
else output << "insert " << name() << ":RhoPhase " << ix << " "
<< _rhophase[ix] << "\n";
}
for(unsigned int ix=0;ix<_rhomass.size();++ix) {
if(ix<1) output << "newdef " << name() << ":RhoMasses " << ix << " "
<< _rhomass[ix]/GeV << "\n";
else output << "insert " << name() << ":RhoMasses " << ix << " "
<< _rhomass[ix]/GeV << "\n";
}
for(unsigned int ix=0;ix<_rhowidth.size();++ix) {
if(ix<1) output << "newdef " << name() << ":RhoWidths " << ix << " "
<< _rhowidth[ix]/GeV << "\n";
else output << "insert " << name() << ":RhoWidths " << ix << " "
<< _rhowidth[ix]/GeV << "\n";
}
// integration weights for the different channels
for(unsigned int ix=0;ix<_zerowgts.size();++ix) {
output << "newdef " << name() << ":AllNeutralWeights "
<< ix << " " << _zerowgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_onewgts.size();++ix) {
output << "newdef " << name() << ":OneChargedWeights "
<< ix << " " << _onewgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_twowgts.size();++ix) {
output << "newdef " << name() << ":TwoChargedWeights "
<< ix << " " << _twowgts[ix] << "\n";
}
for(unsigned int ix=0;ix<_threewgts.size();++ix) {
output << "newdef " << name() << ":ThreeChargedWeights "
<< ix << " " << _threewgts[ix] << "\n";
}
output << "newdef " << name() << ":ZeroMax " << _zeromax << "\n";
output << "newdef " << name() << ":OneMax " << _onemax << "\n";
output << "newdef " << name() << ":TwoMax " << _twomax << "\n";
output << "newdef " << name() << ":ThreeMax " << _threemax << "\n";
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:26 PM (1 d, 10 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805823
Default Alt Text
(67 KB)

Event Timeline