Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F8723782
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
69 KB
Subscribers
None
View Options
diff --git a/Decay/VectorMeson/a1ThreePionCLEODecayer.cc b/Decay/VectorMeson/a1ThreePionCLEODecayer.cc
--- a/Decay/VectorMeson/a1ThreePionCLEODecayer.cc
+++ b/Decay/VectorMeson/a1ThreePionCLEODecayer.cc
@@ -1,1126 +1,1074 @@
// -*- C++ -*-
//
// a1ThreePionCLEODecayer.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 a1ThreePionCLEODecayer class.
//
#include "a1ThreePionCLEODecayer.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/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
using Constants::pi;
a1ThreePionCLEODecayer::a1ThreePionCLEODecayer()
: _rhomass(2), _rhowidth(2), _f2mass(1.275*GeV), _f2width(0.185*GeV),
_pf2cc(ZERO), _pf200(ZERO), _f0mass(1.186*GeV), _f0width(0.350*GeV),
_pf0cc(ZERO), _pf000(ZERO), _sigmamass(0.860*GeV), _sigmawidth(0.880*GeV),
_psigmacc(ZERO), _psigma00(ZERO), _mpi0(ZERO), _mpic(ZERO),
_coupling(45.57/GeV), _rhomagP(2), _rhophaseP(2), _rhomagD(2),
_rhophaseD(2), _f2mag(0.71/GeV2), _f2phase(0.56*pi), _f2coup(ZERO,ZERO),
_f0mag(0.77), _f0phase(-0.54*pi), _f0coup(0.,0.), _sigmamag(2.10),
_sigmaphase(0.23*pi), _sigmacoup(0.,0.), _localparameters(true),
_zerowgts(9), _onewgts(9), _twowgts(9), _threewgts(12), _zeromax(13.0704),
_onemax(6.91104), _twomax(6.94654), _threemax(6.40086) {
// rho masses and widths
_rhomass[0] = 0.7743*GeV; _rhowidth[0] = 0.1491*GeV;
_rhomass[1] = 1.370 *GeV; _rhowidth[1] = 0.386 *GeV;
// p-wave rho and rho prime
_rhomagP[0] = 1. ; _rhophaseP[0] = 0.;
_rhomagP[1] = 0.12; _rhophaseP[1] = 0.99*pi;
// d-wave rho and rho prime
_rhomagD[0] = 0.37/GeV2; _rhophaseD[0] = -0.15*pi;
_rhomagD[1] = 0.87/GeV2; _rhophaseD[1] = 0.53*pi;
// set up the integration channels
_zerowgts[0] = 0.132162;_zerowgts[1] = 0.116638;_zerowgts[2] = 0.121088;
_zerowgts[3] = 0.10656 ;_zerowgts[4] = 0.102577;_zerowgts[5] = 0.101169;
_zerowgts[6] = 0.104587;_zerowgts[7] = 0.104663;_zerowgts[8] = 0.110557;
_onewgts[0] = 0.177017;_onewgts[1] = 0.176011;_onewgts[2] = 0.110129;
_onewgts[3] = 0.108023;_onewgts[4] = 0.110553;_onewgts[5] = 0.109976;
_onewgts[6] = 0.088634;_onewgts[7] = 0.059104;_onewgts[8] = 0.060553;
_twowgts[0] = 0.173357;_twowgts[1] = 0.172283;_twowgts[2] = 0.116031;
_twowgts[3] = 0.114642;_twowgts[4] = 0.109058;_twowgts[5] = 0.114073;
_twowgts[6] = 0.080946;_twowgts[7] = 0.060135;_twowgts[8] = 0.059477;
_threewgts[0] = 0.125022;_threewgts[1] = 0.129911;_threewgts[2] = 0.074165;
_threewgts[3] = 0.075813;_threewgts[4 ]= 0.071154;_threewgts[5 ]= 0.077730;
_threewgts[6] = 0.082255;_threewgts[7 ]= 0.086761;_threewgts[8 ]= 0.067106;
_threewgts[9] = 0.070171;_threewgts[10]= 0.070146;_threewgts[11]= 0.069767;
// generation of intermediates
generateIntermediates(true);
}
void a1ThreePionCLEODecayer::doinit() {
- DecayIntegrator::doinit();
+ DecayIntegrator2::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);
// 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);
// the f_2
tPDPtr f2=getParticleData(225);
// the f_0
tPDPtr f0=getParticleData(10221);
// set up the integration channels
- tPDVector extpart(4);
- DecayPhaseSpaceChannelPtr newchannel;
- DecayPhaseSpaceModePtr mode;
// decay mode a_10 -> pi0 pi0 pi0
- extpart[0]=a10;
- extpart[1]=pi0;
- extpart[2]=pi0;
- extpart[3]=pi0;
- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
+ tPDPtr in = a10;
+ tPDVector out = {pi0,pi0,pi0};
+ PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_zeromax));
+ vector<PhaseSpaceChannel> channels;
// there are six sigma channels
- tPDPtr temp;
for(unsigned int ix=0;ix<3;++ix) {
+ tPDPtr temp;
if(ix==0) temp = sigma;
else if(ix==1) temp = f2;
else if(ix==2) temp = f0;
if(temp) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,1);
- newchannel->addIntermediate(temp,0,0.0,2,3);
- mode->addChannel(newchannel);
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,2);
- newchannel->addIntermediate(temp,0,0.0,1,3);
- mode->addChannel(newchannel);
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,3);
- newchannel->addIntermediate(temp,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,temp,0,1,1,2,1,3));
+ channels.push_back((PhaseSpaceChannel(mode),0,temp,0,2,1,1,1,3));
+ channels.push_back((PhaseSpaceChannel(mode),0,temp,0,3,1,1,1,2));
}
}
- if(_zerowgts.size()!=mode->numberChannels())
- _zerowgts=vector<double>(mode->numberChannels(),1./mode->numberChannels());
- addMode(mode,_zeromax,_zerowgts);
+ if(_zerowgts.size()!=channels.size())
+ _zerowgts=vector<double>(channels.size(),
+ 1./channels.size());
+ for(unsigned int ix=0;ix<channels.size();++ix) {
+ channels[ix].weight(_zerowgts[ix]);
+ mode->addChannel(channels[ix]);
+ }
+ addMode(mode);
// decay mode a_1+ -> pi+ pi0 pi0
- extpart[0]=a1p;
- extpart[1]=pi0;
- extpart[2]=pi0;
- extpart[3]=pip;
- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
+ in = a1p;
+ out = {pi0,pi0,pip};
+ mode = new_ptr(PhaseSpaceMode(in,out,_onemax));
+ channels.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho+ channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,1);
- newchannel->addIntermediate(rhop[ix],0,0.0,2,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rhop[ix],0,1,1,2,1,3));
// second rho+ channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,2);
- newchannel->addIntermediate(rhop[ix],0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rhop[ix],0,2,1,1,1,3));
}
// the sigma channel
if(sigma) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,3);
- newchannel->addIntermediate(sigma,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,sigma,0,3,1,1,1,2));
}
// the f_2 channel
if(f2) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,3);
- newchannel->addIntermediate(f2,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f2,0,3,1,1,1,2));
}
// the f_0 channel
if(f0) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,3);
- newchannel->addIntermediate(f0,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f0,0,3,1,1,1,2));
}
- if(_onewgts.size()!=mode->numberChannels())
- _onewgts=vector<double>(mode->numberChannels(),1./mode->numberChannels());
- addMode(mode,_onemax,_onewgts);
+ if(_onewgts.size()!=channels.size())
+ _onewgts=vector<double>(channels.size(),1./channels.size());
+ for(unsigned int ix=0;ix<channels.size();++ix) {
+ channels[ix].weight(_onewgts[ix]);
+ mode->addChannel(channels[ix]);
+ }
+ addMode(mode);
// decay mode a_10 -> pi+ pi- pi0
- extpart[0]=a10;
- extpart[1]=pip;
- extpart[2]=pim;
- extpart[3]=pi0;
- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
+ in = a10;
+ out = {pip,pim,pi0};
+ mode = new_ptr(PhaseSpaceMode(in,out,_twomax));
+ channels.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(!rhop[ix]) continue;
// first rho channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,1);
- newchannel->addIntermediate(rhom[ix],0,0.0,2,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rhom[ix],0,1,1,2,1,3));
// second channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,2);
- newchannel->addIntermediate(rhop[ix],0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rhom[ix],0,2,1,1,1,3));
}
// sigma channel
if(sigma) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,3);
- newchannel->addIntermediate(sigma,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,sigma,0,3,1,1,1,2));
}
// f_2 channel
if(f2) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,3);
- newchannel->addIntermediate(f2,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f2,0,3,1,1,1,2));
}
// f_0 channel
if(f0) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a10,0,0.0,-1,3);
- newchannel->addIntermediate(f0,0,0.0,1,2);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f0,0,3,1,1,1,2));
}
- if(_twowgts.size()!=mode->numberChannels())
- _twowgts=vector<double>(mode->numberChannels(),1./mode->numberChannels());
- addMode(mode,_twomax,_twowgts);
+ if(_twowgts.size()!=channels.size())
+ _twowgts=vector<double>(channels.size(),1./channels.size());
+ for(unsigned int ix=0;ix<channels.size();++ix) {
+ channels[ix].weight(_twowgts[ix]);
+ mode->addChannel(channels[ix]);
+ }
+ addMode(mode);
// decay mode a_1+ -> pi+ pi+ pi-
- extpart[0]=a1p;
- extpart[1]=pip;
- extpart[2]=pip;
- extpart[3]=pim;
- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
+ in = a1p;
+ out = {pip,pip,pim};
+ mode = new_ptr(PhaseSpaceMode(in,out,_threemax));
+ channels.clear();
for(unsigned int ix=0;ix<3;++ix) {
// the neutral rho channels
if(!rho0[ix]) continue;
// first channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,1);
- newchannel->addIntermediate(rho0[ix],0,0.0,2,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rho0[ix],0,1,1,2,1,3));
// interchanged channel
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,2);
- newchannel->addIntermediate(rho0[ix],0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,rho0[ix],0,2,1,1,1,3));
}
// the sigma channels
if(sigma) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,1);
- newchannel->addIntermediate(sigma,0,0.0,2,3);
- mode->addChannel(newchannel);
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,2);
- newchannel->addIntermediate(sigma,0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,sigma,0,1,1,2,1,3));
+ channels.push_back((PhaseSpaceChannel(mode),0,sigma,0,2,1,1,1,3));
}
// the f_2 channels
if(f2) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,1);
- newchannel->addIntermediate(f2,0,0.0,2,3);
- mode->addChannel(newchannel);
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,2);
- newchannel->addIntermediate(f2,0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f2,0,1,1,2,1,3));
+ channels.push_back((PhaseSpaceChannel(mode),0,f2,0,2,1,1,1,3));
}
// the f_0 channel
if(f0) {
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,1);
- newchannel->addIntermediate(f0,0,0.0,2,3);
- mode->addChannel(newchannel);
- newchannel = new_ptr(DecayPhaseSpaceChannel(mode));
- newchannel->addIntermediate(a1p,0,0.0,-1,2);
- newchannel->addIntermediate(f0,0,0.0,1,3);
- mode->addChannel(newchannel);
+ channels.push_back((PhaseSpaceChannel(mode),0,f0,0,1,1,2,1,3));
+ channels.push_back((PhaseSpaceChannel(mode),0,f0,0,2,1,1,1,3));
}
- if(_threewgts.size()!=mode->numberChannels())
- _threewgts=vector<double>(mode->numberChannels(),1./mode->numberChannels());
- addMode(mode,_threemax,_threewgts);
+ if(_threewgts.size()!=channels.size())
+ _threewgts=vector<double>(channels.size(),1./channels.size());
+ for(unsigned int ix=0;ix<channels.size();++ix) {
+ channels[ix].weight(_threewgts[ix]);
+ mode->addChannel(channels[ix]);
+ }
+ 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]);
}
resetIntermediate(sigma,_sigmamass,_sigmawidth);
resetIntermediate(f2,_f2mass,_f2width);
resetIntermediate(f0,_f0mass,_f0width);
// 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) {
_rhomass[ix]=rhop[ix]->mass();
_rhowidth[ix]=rhop[ix]->width();
}
if(f2) {
_f2mass=f2->mass();
_f2width=f2->width();
}
if(f0) {
_f0mass=f0->mass();
_f0width=f0->width();
}
if(sigma) {
_sigmamass=sigma->mass();
_sigmawidth=sigma->width();
}
}
// parameters for the breit-wigners
_mpic=pip->mass();
_mpi0=pi0->mass();
// momenta of the decay products for on-shell particles
_psigmacc = Kinematics::pstarTwoBodyDecay(_sigmamass,_mpic,_mpic);
_psigma00 = Kinematics::pstarTwoBodyDecay(_sigmamass,_mpi0,_mpi0);
_pf2cc = Kinematics::pstarTwoBodyDecay(_f2mass ,_mpic,_mpic);
_pf200 = Kinematics::pstarTwoBodyDecay(_f2mass ,_mpi0,_mpi0);
_pf0cc = Kinematics::pstarTwoBodyDecay(_f0mass ,_mpic,_mpic);
_pf000 = Kinematics::pstarTwoBodyDecay(_f0mass ,_mpi0,_mpi0);
_prhocc.resize(3);_prhoc0.resize(3);
for(unsigned int ix=0;ix<3;++ix) {
_prhocc[ix] = Kinematics::pstarTwoBodyDecay(_rhomass[ix],_mpic,_mpic);
_prhoc0[ix] = Kinematics::pstarTwoBodyDecay(_rhomass[ix],_mpic,_mpi0);
}
// couplings for the different modes
Complex ii(0.,1.);
_rhocoupP.resize(_rhomagP.size());
for(unsigned int ix=0;ix<_rhomagP.size();++ix)
_rhocoupP[ix]=_rhomagP[ix]*(cos(_rhophaseP[ix])+ii*sin(_rhophaseP[ix]));
_rhocoupD.resize(_rhomagD.size());
for(unsigned int ix=0;ix<_rhomagD.size();++ix)
_rhocoupD[ix]=_rhomagD[ix]*(cos(_rhophaseD[ix])+ii*sin(_rhophaseD[ix]));
_f0coup=_f0mag*(cos(_f0phase)+ii*sin(_f0phase));
_f2coup=_f2mag*(cos(_f2phase)+ii*sin(_f2phase));
_sigmacoup=_sigmamag*(cos(_sigmaphase)+ii*sin(_sigmaphase));
}
inline void a1ThreePionCLEODecayer::doinitrun() {
- DecayIntegrator::doinitrun();
+ DecayIntegrator2::doinitrun();
if(initialize()) {
// get the weights for the different channels
for(unsigned int ix=0;ix<_zerowgts.size();++ix)
- _zerowgts[ix]=mode(0)->channelWeight(ix);
+ _zerowgts[ix]=mode(0)->channels()[ix].weight();
for(unsigned int ix=0;ix<_onewgts.size();++ix)
- _onewgts[ix]=mode(1)->channelWeight(ix);
+ _onewgts[ix]=mode(1)->channels()[ix].weight();
for(unsigned int ix=0;ix<_twowgts.size();++ix)
- _twowgts[ix]=mode(2)->channelWeight(ix);
+ _twowgts[ix]=mode(2)->channels()[ix].weight();
for(unsigned int ix=0;ix<_threewgts.size();++ix)
- _threewgts[ix]=mode(3)->channelWeight(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();
}
}
int a1ThreePionCLEODecayer::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 a1ThreePionCLEODecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_rhomass,GeV) << ounit(_rhowidth,GeV)
<< ounit(_prhocc,GeV) << ounit(_prhoc0,GeV)
<< ounit(_f2mass,GeV) << ounit(_f2width,GeV)
<< ounit(_pf2cc,GeV)
<< ounit(_pf200,GeV) << ounit(_f0mass,GeV)
<< ounit(_f0width,GeV) << ounit(_pf0cc,GeV) << ounit(_pf000,GeV)
<< ounit(_sigmamass,GeV) << ounit(_sigmawidth,GeV)
<< ounit(_psigmacc,GeV) << ounit(_psigma00,GeV)
<< ounit(_mpi0,GeV) << ounit(_mpic,GeV)
<< ounit(_coupling,1/GeV) << _rhomagP << _rhophaseP
<< _rhocoupP << ounit(_rhomagD ,1/GeV2)<< _rhophaseD
<< ounit(_rhocoupD,1/GeV2) << ounit(_f2mag,1/GeV2)
<< _f2phase << ounit(_f2coup,1/GeV2)
<< _f0mag << _f0phase << _f0coup
<< _sigmamag << _sigmaphase << _sigmacoup
<< _localparameters
<< _zerowgts << _onewgts << _twowgts << _threewgts
<< _zeromax << _onemax << _twomax << _threemax;
}
void a1ThreePionCLEODecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_rhomass,GeV) >> iunit(_rhowidth,GeV)
>> iunit(_prhocc,GeV) >> iunit(_prhoc0,GeV)
>> iunit(_f2mass,GeV) >> iunit(_f2width,GeV)
>> iunit(_pf2cc,GeV)
>> iunit(_pf200,GeV) >> iunit(_f0mass,GeV)
>> iunit(_f0width,GeV) >> iunit(_pf0cc,GeV) >> iunit(_pf000,GeV)
>> iunit(_sigmamass,GeV) >> iunit(_sigmawidth,GeV)
>> iunit(_psigmacc,GeV) >> iunit(_psigma00,GeV)
>> iunit(_mpi0,GeV) >> iunit(_mpic,GeV)
>> iunit(_coupling,1/GeV) >> _rhomagP >> _rhophaseP
>> _rhocoupP >> iunit(_rhomagD,1/GeV2) >> _rhophaseD
>> iunit(_rhocoupD,1/GeV2)>>iunit(_f2mag,1/GeV2)
>> _f2phase >> iunit(_f2coup,1/GeV2)
>> _f0mag >> _f0phase >> _f0coup
>> _sigmamag >> _sigmaphase >> _sigmacoup
>> _localparameters
>> _zerowgts >> _onewgts >> _twowgts >> _threewgts
>> _zeromax >> _onemax >> _twomax >> _threemax;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<a1ThreePionCLEODecayer,DecayIntegrator>
+DescribeClass<a1ThreePionCLEODecayer,DecayIntegrator2>
describeHerwiga1ThreePionCLEODecayer("Herwig::a1ThreePionCLEODecayer", "HwVMDecay.so");
void a1ThreePionCLEODecayer::Init() {
static ClassDocumentation<a1ThreePionCLEODecayer> documentation
("The a1ThreePionCLEODecayer class performs the decay of the "
"a_1 to three pions using the model of CLEO",
"The decay of a_1 to three pions was modelled after \\cite{Asner:1999kj}.",
"%\\cite{Asner:1999kj}\n"
"\\bibitem{Asner:1999kj}\n"
" D.~M.~Asner {\\it et al.} [CLEO Collaboration],\n"
" ``Hadronic structure in the decay tau- --> nu/tau pi- pi0 pi0 and the sign\n"
" %of the tau neutrino helicity,''\n"
" Phys.\\ Rev.\\ D {\\bf 61}, 012002 (2000)\n"
" [arXiv:hep-ex/9902022].\n"
" %%CITATION = PHRVA,D61,012002;%%\n"
);
static ParVector<a1ThreePionCLEODecayer,Energy> interfacerhomass
("RhoMasses",
"The masses of the different rho resonnaces",
&a1ThreePionCLEODecayer::_rhomass,
GeV, 0, ZERO, -10000*GeV, 10000*GeV, false, false, true);
static ParVector<a1ThreePionCLEODecayer,Energy> interfacerhowidth
("RhoWidths",
"The widths of the different rho resonnaces",
&a1ThreePionCLEODecayer::_rhowidth,
GeV, 0, ZERO, -10000*GeV, 10000*GeV, false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacef_2Mass
("f_2Mass",
"The mass of the f_2 meson",
&a1ThreePionCLEODecayer::_f2mass, GeV, 1.275*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacef_2Width
("f_2Width",
"The width of the f_2 meson",
&a1ThreePionCLEODecayer::_f2width, GeV, 0.185*GeV, ZERO, 1.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacef_0Mass
("f_0Mass",
"The mass of the f_0 meson",
&a1ThreePionCLEODecayer::_f0mass, GeV, 1.186*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacef_0Width
("f_0Width",
"The width of the f_0 meson",
&a1ThreePionCLEODecayer::_f0width, GeV, 0.350*GeV, ZERO, 1.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacesigmaMass
("sigmaMass",
"The mass of the sigma meson",
&a1ThreePionCLEODecayer::_sigmamass, GeV, 0.860*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,Energy> interfacesigmaWidth
("sigmaWidth",
"The width of the sigma meson",
&a1ThreePionCLEODecayer::_sigmawidth, GeV, 0.880*GeV, ZERO, 2.0*GeV,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,InvEnergy> interfaceCoupling
("Coupling",
"The overall coupling for the decay",
&a1ThreePionCLEODecayer::_coupling, 1./GeV, 45.57/GeV, -0./GeV, 1000./GeV,
false, false, false);
static ParVector<a1ThreePionCLEODecayer,double> interfacerhomagP
("RhoPWaveMagnitude",
"The magnitude of the couplings for the p-wave rho currents",
&a1ThreePionCLEODecayer::_rhomagP,
0, 0, 0, 0, 10000, false, false, true);
static ParVector<a1ThreePionCLEODecayer,double> interfacerhophaseP
("RhoPWavePhase",
"The phase of the couplings for the p-wave rho currents",
&a1ThreePionCLEODecayer::_rhophaseP,
0, 0, 0, -pi, pi, false, false, true);
static ParVector<a1ThreePionCLEODecayer,InvEnergy2> interfacerhomagD
("RhoDWaveMagnitude",
"The magnitude of the couplings for the d-wave rho currents",
&a1ThreePionCLEODecayer::_rhomagD,
1/MeV2, 0, ZERO, ZERO, 10000/MeV2, false, false, true);
static ParVector<a1ThreePionCLEODecayer,double> interfacerhophaseD
("RhoDWavePhase",
"The phase of the couplings for the d-wave rho currents",
&a1ThreePionCLEODecayer::_rhophaseD,
0, 0, 0, -pi, pi, false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfacef0Phase
("f0Phase",
"The phase of the f_0 scalar current",
&a1ThreePionCLEODecayer::_f0phase, 0.54*pi, -pi, pi,
false, false, Interface::limited);
static Parameter<a1ThreePionCLEODecayer,double> interfacef2Phase
("f2Phase",
"The phase of the f_2 tensor current",
&a1ThreePionCLEODecayer::_f2phase, 0.56*pi, -pi, pi,
false, false, Interface::limited);
static Parameter<a1ThreePionCLEODecayer,double> interfacesigmaPhase
("sigmaPhase",
"The phase of the sigma scalar current",
&a1ThreePionCLEODecayer::_sigmaphase, 0.23*pi, -pi, pi,
false, false, Interface::limited);
static Parameter<a1ThreePionCLEODecayer,double> interfacef0Magnitude
("f0Magnitude",
"The magnitude of the f_0 scalar current",
&a1ThreePionCLEODecayer::_f0mag, 0.77, 0.0, 10,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,InvEnergy2> interfacef2Magnitude
("f2Magnitude",
"The magnitude of the f_2 tensor current",
&a1ThreePionCLEODecayer::_f2mag, 1./GeV2, 0.71/GeV2, ZERO, 10./GeV2,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfacesigmaMagnitude
("sigmaMagnitude",
"The magnitude of the sigma scalar current",
&a1ThreePionCLEODecayer::_sigmamag, 2.1, 0.0, 10,
false, false, true);
static Switch<a1ThreePionCLEODecayer,bool> interfaceLocalParameters
("LocalParameters",
"Use local values of the intermediate resonances masses and widths",
&a1ThreePionCLEODecayer::_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 ParVector<a1ThreePionCLEODecayer,double> interfacezerowgts
("AllNeutralWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^0->pi0pi0pi0",
&a1ThreePionCLEODecayer::_zerowgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionCLEODecayer,double> interfaceonewgts
("OneChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi0pi0",
&a1ThreePionCLEODecayer::_onewgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionCLEODecayer,double> interfacetwowgts
("TwoChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^0->pi+pi-pi0",
&a1ThreePionCLEODecayer::_twowgts,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<a1ThreePionCLEODecayer,double> interfacethreewgts
("ThreeChargedWeights",
"The weights of the different channels to use for the integration of"
" the decay a_1^+->pi+pi+pi-",
&a1ThreePionCLEODecayer::_threewgts,
0, 0, 0, -10000, 10000, false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfaceZeroMax
("ZeroMax",
"The maximum weight for the integration fo the channel a_1^0->pi0pi0pi0",
&a1ThreePionCLEODecayer::_zeromax, 0.0716349E3, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfaceOneMax
("OneMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi0pi0",
&a1ThreePionCLEODecayer::_onemax,1.23756E3, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfaceTwoMax
("TwoMax",
"The maximum weight for the integration fo the channel a_1^0->pi+pi-pi0",
&a1ThreePionCLEODecayer::_twomax,2.43819E3, 0.0, 10000.0,
false, false, true);
static Parameter<a1ThreePionCLEODecayer,double> interfaceThreeMax
("ThreeMax",
"The maximum weight for the integration fo the channel a_1^+->pi+pi+pi-",
&a1ThreePionCLEODecayer::_threemax, 1.38754E3, 0.0, 10000.0,
false, false, true);
}
-double a1ThreePionCLEODecayer::me2(const int ichan,
- const Particle & inpart,
- const ParticleVector & decay,
- MEOption meopt) const {
+void a1ThreePionCLEODecayer::
+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 a1ThreePionCLEODecayer::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>(&inpart),
+ const_ptr_cast<tPPtr>(&part),
incoming,false);
}
- if(meopt==Terminate) {
- VectorWaveFunction::constructSpinInfo(_vectors,const_ptr_cast<tPPtr>(&inpart),
- 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);
- return 0.;
- }
// momentum of the incoming particle
- Lorentz5Momentum Q=inpart.momentum();
+ Lorentz5Momentum Q=part.momentum();
Energy2 q2=Q.mass2();
// identify the mesons
// calculate the invariants and form factors
- Lorentz5Momentum ps1=decay[1]->momentum()+decay[2]->momentum();
- Lorentz5Momentum ps2=decay[0]->momentum()+decay[2]->momentum();
- Lorentz5Momentum ps3=decay[0]->momentum()+decay[1]->momentum();
+ Lorentz5Momentum ps1=momenta[1]+momenta[2];
+ Lorentz5Momentum ps2=momenta[0]+momenta[2];
+ Lorentz5Momentum ps3=momenta[0]+momenta[1];
ps1.rescaleMass();
ps2.rescaleMass();
ps3.rescaleMass();
Energy2 s1=ps1.mass2(),s2=ps2.mass2(),s3=ps3.mass2();
complex<InvEnergy> F1,F2,F3;
formFactors(imode(),ichan,q2,s1,s2,s3,F1,F2,F3);
// use the form-factors to compute the current
LorentzPolarizationVector output=
- F1*decay[1]->momentum()-F1*decay[2]->momentum()
- -F2*decay[2]->momentum()+F2*decay[0]->momentum()
- +F3*decay[0]->momentum()-F3*decay[1]->momentum();
+ F1*momenta[1]-F1*momenta[2]
+ -F2*momenta[2]+F2*momenta[0]
+ +F3*momenta[0]-F3*momenta[1];
// compute the matrix element
for(unsigned int ix=0;ix<3;++ix)
(*ME())(ix,0,0,0)=output.dot(_vectors[ix]);
// answer
double out = ME()->contract(_rho).real();
// test of the answer
-// double test = threeBodyMatrixElement(imode(),sqr(inpart.mass()),
-// s3,s2,s1,decay[0]->mass(),decay[1]->mass(),
-// decay[2]->mass());
-// if(ichan<0) cerr << "testing matrix element " << inpart.PDGName() << " -> "
-// << decay[0]->PDGName() << " " << decay[1]->PDGName() << " "
-// << decay[2]->PDGName() << out << " " << test << " "
+// 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 the answer
return out;
}
// matrix element for the running a_1 width
double a1ThreePionCLEODecayer::
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 m12=m1*m1,m22=m2*m2,m32=m3*m3;
// calculate the form factors
complex<InvEnergy> F1,F2,F3;
formFactors(iopt,-1,q2,s1,s2,s3,F1,F2,F3);
// analytic calculation of the matrix element
double dot1=( F1*conj(F1)*(2.*m22+2.*m32-s1)+F2*conj(F2)*(2.*m12+2.*m32-s2)
+F3*conj(F3)*(2.*m12+2.*m22-s3)-F1*conj(F2)*( s1+s2-s3-4.*m32)
+F1*conj(F3)*( s1-s2+s3-4.*m22)-F2*conj(F3)*(-s1+s2+s3-4.*m12)).real();
complex<Energy> dot2 = 0.5*(F1*(s3-m32-s2+m22)-F2*(s1-m12-s3+m32)+F3*(s2-m22-s1+m12));
return (-dot1+(dot2*conj(dot2)).real()/q2)/3.;
}
WidthCalculatorBasePtr
a1ThreePionCLEODecayer::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;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 m[3];
m[0] = ncharged<2 ? _mpi0 : _mpic;
m[1] = m[0];
m[2] = (ncharged==0||ncharged==2) ? _mpi0 : _mpic;
return new_ptr(ThreeBodyAllOnCalculator<a1ThreePionCLEODecayer>
(inweights,intype,inmass,inwidth,inpow,*this,ncharged,m[0],m[1],m[2]));
}
// calculate the form factos
void a1ThreePionCLEODecayer::formFactors(int iopt,int ichan,
Energy2 q2,Energy2 s1,Energy2 s2,
Energy2 s3,
complex<InvEnergy> & FF1,
complex<InvEnergy> & FF2,
complex<InvEnergy> & FF3) const {
Complex F1, F2, F3;
InvEnergy fact = _coupling;
// a_1^0 pi0 pi0 pi0 mode
if(iopt==0) {
fact*=1./sqrt(6.);
// compute the breit wigners we need
Complex sigbws1 = sigmaBreitWigner(s1,1);
Complex sigbws2 = sigmaBreitWigner(s2,1);
Complex sigbws3 = sigmaBreitWigner(s3,1);
Complex f0bws1 = f0BreitWigner(s1,1);
Complex f0bws2 = f0BreitWigner(s2,1);
Complex f0bws3 = f0BreitWigner(s3,1);
Complex f2bws1 = f2BreitWigner(s1,1);
Complex f2bws2 = f2BreitWigner(s2,1);
Complex f2bws3 = f2BreitWigner(s3,1);
if(ichan<0) {
// the scalar terms
F1=2./3.*(_sigmacoup*sigbws3+_f0coup*f0bws3)
-2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2);
F2=2./3.*(_sigmacoup*sigbws3+_f0coup*f0bws3)
-2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1);
F3=-2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1)
+2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2);
// the tensor terms
complex<Energy2> Dfact1 = 1./18.*(4.*_mpi0*_mpi0-s1)*(q2+s1-_mpi0*_mpi0)/s1*f2bws1;
complex<Energy2> Dfact2 = 1./18.*(4.*_mpi0*_mpi0-s2)*(q2+s2-_mpi0*_mpi0)/s2*f2bws2;
complex<Energy2> Dfact3 = 1./18.*(4.*_mpi0*_mpi0-s3)*(q2-_mpi0*_mpi0+s3)/s3*f2bws3;
F1+=_f2coup*( 0.5*(s3-s2)*f2bws1-Dfact2+Dfact3);
F2+=_f2coup*( 0.5*(s3-s1)*f2bws2-Dfact1+Dfact3);
F3+=_f2coup*(-0.5*(s1-s2)*f2bws3-Dfact1+Dfact2);
}
else if(ichan==0) {
F2=-2./3.*_sigmacoup*sigbws1;
F3=-2./3.*_sigmacoup*sigbws1;
}
else if(ichan==1) {
F1=-2./3.*_sigmacoup*sigbws2;
F3=+2./3.*_sigmacoup*sigbws2;
}
else if(ichan==2) {
F1= 2./3.*_sigmacoup*sigbws3;
F2= 2./3.*_sigmacoup*sigbws3;
}
else if(ichan==3) {
complex<Energy2> Dfact1 = 1./18.*(4.*_mpi0*_mpi0-s1)*(q2+s1-_mpi0*_mpi0)/s1*f2bws1;
F1+=_f2coup*0.5*(s3-s2)*f2bws1;
F2-=_f2coup*Dfact1;
F3-=_f2coup*Dfact1;
}
else if(ichan==4) {
complex<Energy2> Dfact2 = 1./18.*(4.*_mpi0*_mpi0-s2)*(q2+s2-_mpi0*_mpi0)/s2*f2bws2;
F2+=_f2coup*0.5*(s3-s1)*f2bws2;
F1-=_f2coup*Dfact2;
F3+=_f2coup*Dfact2;
}
else if(ichan==5) {
complex<Energy2> Dfact3 = 1./18.*(4.*_mpi0*_mpi0-s3)*(q2-_mpi0*_mpi0+s3)/s3*f2bws3;
F3+=-_f2coup*0.5*(s1-s2)*f2bws3;
F1+=_f2coup*Dfact3;
F2+=_f2coup*Dfact3;
}
else if(ichan==6) {
F2=-2./3.*_f0coup*f0bws1;
F3=-2./3.*_f0coup*f0bws1;
}
else if(ichan==7) {
F1=-2./3.*_f0coup*f0bws2;
F3=+2./3.*_f0coup*f0bws2;
}
else if(ichan==8) {
F1= 2./3.*_f0coup*f0bws3;
F2= 2./3.*_f0coup*f0bws3;
}
}
// a_1^+ -> pi0 pi0 pi+
else if(iopt==1) {
fact *= 1./sqrt(2.);
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3],f0bw,sigbw,f2bw;
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix]=rhoBreitWigner(ix,s1,1);
rhos2bw[ix]=rhoBreitWigner(ix,s2,1);
}
f0bw = f0BreitWigner(s3,1);
sigbw = sigmaBreitWigner(s3,1);
f2bw = f2BreitWigner(s3,1);
if(ichan<0) {
// the p-wave rho terms
for(unsigned int ix=0;ix<_rhocoupP.size();++ix) {
F1+=_rhocoupP[ix]*rhos1bw[ix];
F2+=_rhocoupP[ix]*rhos2bw[ix];
}
// the D-wave rho terms
Energy2 Dfact1=-1./3.*((s3-_mpic*_mpic)-(s1-_mpi0*_mpi0));
Energy2 Dfact2=-1./3.*((s3-_mpic*_mpic)-(s2-_mpi0*_mpi0));
for(unsigned int ix=0;ix<_rhocoupD.size();++ix) {
F1+=Dfact1*_rhocoupD[ix]*rhos2bw[ix];
F2+=Dfact2*_rhocoupD[ix]*rhos1bw[ix];
F3+=_rhocoupD[ix]*(Dfact2*rhos1bw[ix]-Dfact1*rhos2bw[ix]);
}
// the scalar terms
Complex scalar=2./3.*(_sigmacoup*sigbw+_f0coup*f0bw);
F1+=scalar;
F2+=scalar;
// the tensor terms
Complex Dfact3=1./18./s3*_f2coup*(q2-_mpic*_mpic+s3)*(4.*_mpi0*_mpi0-s3)*f2bw;
F1+=Dfact3;F2+=Dfact3;
F3-=0.5*_f2coup*(s1-s2)*f2bw;
}
else if(ichan%2==0&&ichan<=4) {
unsigned int ires=ichan/2;
if(ires<_rhocoupP.size()) F1+=_rhocoupP[ires]*rhos1bw[ires];
Energy2 Dfact2=-1./3.*((s3-_mpic*_mpic)-(s2-_mpi0*_mpi0));
if(ires<_rhocoupD.size()) {
F2+=Dfact2*_rhocoupD[ires]*rhos1bw[ires];
F3+=_rhocoupD[ires]*Dfact2*rhos1bw[ires];
}
}
else if(ichan%2==1&&ichan<=5) {
unsigned int ires=(ichan-1)/2;
if(ires<_rhocoupP.size()) F2+=_rhocoupP[ires]*rhos2bw[ires];
Energy2 Dfact1=-1./3.*((s3-_mpic*_mpic)-(s1-_mpi0*_mpi0));
if(ires<_rhocoupD.size()) {
F1+=Dfact1*_rhocoupD[ires]*rhos2bw[ires];
F3-=_rhocoupD[ires]*Dfact1*rhos2bw[ires];
}
}
else if(ichan==6) {
F1+=2./3.*_sigmacoup*sigbw;
F2+=2./3.*_sigmacoup*sigbw;
}
else if(ichan==7) {
Complex Dfact3=1./18./s3*_f2coup*(q2-_mpic*_mpic+s3)*(4.*_mpi0*_mpi0-s3)*f2bw;
F1+=Dfact3;
F2+=Dfact3;
F3-=0.5*_f2coup*(s1-s2)*f2bw;
}
else if(ichan==8) {
F1+=2./3.*_f0coup*f0bw;
F2+=2./3.*_f0coup*f0bw;
}
}
// a_1^0 ->pi+pi-pi0
else if(iopt==2) {
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3],f0bw,sigbw,f2bw;
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix]=rhoBreitWigner(ix,s1,1);
rhos2bw[ix]=rhoBreitWigner(ix,s2,1);
}
f0bw =f0BreitWigner(s3,0);
sigbw =sigmaBreitWigner(s3,0);
f2bw =f2BreitWigner(s3,0);
if(ichan<0) {
// the p-wave rho terms
for(unsigned int ix=0;ix<_rhocoupP.size();++ix) {
F1+=_rhocoupP[ix]*rhos1bw[ix];
F2+=_rhocoupP[ix]*rhos2bw[ix];
}
// the D-wave rho terms
Energy2 Dfact1=-1./3.*(s3-_mpi0*_mpi0-s1+_mpic*_mpic);
Energy2 Dfact2=-1./3.*(s3-_mpi0*_mpi0-s2+_mpic*_mpic);
for(unsigned int ix=0;ix<_rhocoupD.size();++ix) {
F1+=Dfact1*_rhocoupD[ix]*rhos2bw[ix];
F2+=Dfact2*_rhocoupD[ix]*rhos1bw[ix];
F3+=_rhocoupD[ix]*(Dfact2*rhos1bw[ix]-Dfact1*rhos2bw[ix]);
}
// the scalar terms
Complex scalar=2./3.*(_sigmacoup*sigbw+_f0coup*f0bw);
F1+=scalar;
F2+=scalar;
// the tensor terms
Complex Dfact3=1./18./s3*_f2coup*(q2-_mpi0*_mpi0+s3)*(4.*_mpic*_mpic-s3)*f2bw;
F1+=Dfact3;
F2+=Dfact3;
F3-=0.5*_f2coup*(s1-s2)*f2bw;
}
else if(ichan%2==0&&ichan<=4) {
unsigned int ires=ichan/2;
if(ires<_rhocoupP.size()) F1+=_rhocoupP[ires]*rhos1bw[ires];
Energy2 Dfact2=-1./3.*(s3-_mpi0*_mpi0-s2+_mpic*_mpic);
if(ires<_rhocoupD.size()) {
F2+=Dfact2*_rhocoupD[ires]*rhos1bw[ires];
F3+=_rhocoupD[ires]*Dfact2*rhos1bw[ires];
}
}
else if(ichan%2==1&&ichan<=5) {
unsigned int ires=(ichan-1)/2;
if(ires<_rhocoupP.size()) F2+=_rhocoupP[ires]*rhos2bw[ires];
Energy2 Dfact1=-1./3.*(s3-_mpi0*_mpi0-s1+_mpic*_mpic);
if(ires<_rhocoupD.size()) {
F1+=Dfact1*_rhocoupD[ires]*rhos2bw[ires];
F3-=_rhocoupD[ires]*-Dfact1*rhos2bw[ires];
}
}
else if(ichan==6) {
F1+=2./3.*_sigmacoup*sigbw;
F2+=2./3.*_sigmacoup*sigbw;
}
else if(ichan==7) {
Complex Dfact3=1./18./s3*_f2coup*(q2-_mpi0*_mpi0+s3)*(4.*_mpic*_mpic-s3)*f2bw;
F1+=Dfact3;
F2+=Dfact3;
F3-=0.5*_f2coup*(s1-s2)*f2bw;
}
else if(ichan==8) {
F1+=2./3.*_f0coup*f0bw;
F2+=2./3.*_f0coup*f0bw;
}
}
// a_1^+ -> pi+ pi+ pi- mode
else {
fact *= 1./sqrt(2.);
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3],f0bws1,sigbws1,f2bws1,f0bws2,sigbws2,f2bws2;
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix]=rhoBreitWigner(ix,s1,0);
rhos2bw[ix]=rhoBreitWigner(ix,s2,0);
}
f0bws1 =f0BreitWigner(s1,0);
sigbws1 =sigmaBreitWigner(s1,0);
f2bws1 =f2BreitWigner(s1,0);
f0bws2 =f0BreitWigner(s2,0);
sigbws2 =sigmaBreitWigner(s2,0);
f2bws2 =f2BreitWigner(s2,0);
if(ichan<0) {
// the p-wave rho terms
for(unsigned int ix=0;ix<_rhocoupP.size();++ix) {
F1-=_rhocoupP[ix]*rhos1bw[ix];
F2-=_rhocoupP[ix]*rhos2bw[ix];
}
// the D-wave rho terms
Energy2 Dfact1=1./3.*(s1-s3);
Energy2 Dfact2=1./3.*(s2-s3);
for(unsigned int ix=0;ix<_rhocoupD.size();++ix) {
F1-=Dfact1*_rhocoupD[ix]*rhos2bw[ix];
F2-=Dfact2*_rhocoupD[ix]*rhos1bw[ix];
F3-=_rhocoupD[ix]*(Dfact2*rhos1bw[ix]-Dfact1*rhos2bw[ix]);
}
// the scalar terms
F1-=2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2);
F2-=2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1);
F3+=-2./3.*(_sigmacoup*sigbws1+_f0coup*f0bws1)
+2./3.*(_sigmacoup*sigbws2+_f0coup*f0bws2);
// the tensor terms
complex<Energy2> sfact1
= 1./18.*(4.*_mpic*_mpic-s1)*(q2+s1-_mpic*_mpic)/s1*f2bws1;
complex<Energy2> sfact2
= 1./18.*(4.*_mpic*_mpic-s2)*(q2+s2-_mpic*_mpic)/s2*f2bws2;
F1+=_f2coup*(0.5*(s3-s2)*f2bws1-sfact2);
F2+=_f2coup*(0.5*(s3-s1)*f2bws2-sfact1);
F3+=_f2coup*(-sfact1+sfact2);
}
else if(ichan%2==0&&ichan<=4) {
unsigned int ires=ichan/2;
Energy2 Dfact2=1./3.*(s2-s3);
if(ires<_rhocoupP.size()) F1-=_rhocoupP[ires]*rhos1bw[ires];
if(ires<_rhocoupD.size()){
F2-=Dfact2*_rhocoupD[ires]*rhos1bw[ires];
F3-=_rhocoupD[ires]*Dfact2*rhos1bw[ires];
}
}
else if(ichan%2==1&&ichan<=5) {
unsigned int ires=(ichan-1)/2;
Energy2 Dfact1=1./3.*(s1-s3);
if(ires<_rhocoupP.size()) F2-=_rhocoupP[ires]*rhos2bw[ires];
if(ires<_rhocoupD.size()) {
F1-=Dfact1*_rhocoupD[ires]*rhos2bw[ires];
F3+=_rhocoupD[ires]*Dfact1*rhos2bw[ires];
}
}
else if(ichan==6) {
F2-=2./3.*_sigmacoup*sigbws1;
F3-=2./3.*_sigmacoup*sigbws1;
}
else if(ichan==7) {
F1-=2./3.*_sigmacoup*sigbws2;
F3+=2./3.*_sigmacoup*sigbws2;
}
else if(ichan==8) {
complex<Energy2> sfact1 = 1./18.*(4.*_mpic*_mpic-s1)*(q2+s1-_mpic*_mpic)/s1*f2bws1;
F1+=_f2coup*0.5*(s3-s2)*f2bws1;
F2-=_f2coup*sfact1;
F3-=_f2coup*sfact1;
}
else if(ichan==9) {
complex<Energy2> sfact2 = 1./18.*(4.*_mpic*_mpic-s2)*(q2+s2-_mpic*_mpic)/s2*f2bws2;
F1-=_f2coup*sfact2;
F2+=_f2coup*0.5*(s3-s1)*f2bws2;
F3+=_f2coup*sfact2;
}
else if(ichan==10) {
F2-=2./3.*_f0coup*f0bws1;
F3-=2./3.*_f0coup*f0bws1;
}
else if(ichan==11) {
F1-=2./3.*_f0coup*f0bws2;
F3+=2./3.*_f0coup*f0bws2;
}
}
FF1 = F1 * fact;
FF2 = F2 * fact;
FF3 = F3 * fact;
}
// output the setup information for the particle database
void a1ThreePionCLEODecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
- // parameters for the DecayIntegrator base class
- DecayIntegrator::dataBaseOutput(output,false);
+ // parameters for the DecayIntegrator2 base class
+ DecayIntegrator2::dataBaseOutput(output,false);
// masses and widths of the intermediate particles
output << "newdef " << name() << ":f_2Mass " << _f2mass/GeV << "\n";
output << "newdef " << name() << ":f_2Width " << _f2width/GeV << "\n";
output << "newdef " << name() << ":f_0Mass " << _f0mass/GeV << "\n";
output << "newdef " << name() << ":f_0Width " << _f0width/GeV << "\n";
output << "newdef " << name() << ":sigmaMass " << _sigmamass/GeV << "\n";
output << "newdef " << name() << ":sigmaWidth " << _sigmawidth/GeV << "\n";
for(unsigned int ix=0;ix<_rhomass.size();++ix) {
if(ix<2) 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<2) output << "newdef " << name() << ":RhoWidths " << ix << " "
<< _rhowidth[ix]/GeV << "\n";
else output << "insert " << name() << ":RhoWidths " << ix << " "
<< _rhowidth[ix]/GeV << "\n";
}
// couplings and phases for different channels
output << "newdef " << name() << ":f0Phase " << _f0phase << "\n";
output << "newdef " << name() << ":f2Phase " << _f2phase<< "\n";
output << "newdef " << name() << ":sigmaPhase " << _sigmaphase<< "\n";
output << "newdef " << name() << ":f0Magnitude " << _f0mag<< "\n";
output << "newdef " << name() << ":f2Magnitude " << _f2mag*GeV2 << "\n";
output << "newdef " << name() << ":sigmaMagnitude " << _sigmamag << "\n";
output << "newdef " << name() << ":Coupling " << _coupling*GeV << "\n";
for(unsigned int ix=0;ix<_rhomagP.size();++ix) {
if(ix<2) output << "newdef " << name() << ":RhoPWaveMagnitude " << ix << " "
<< _rhomagP[ix] << "\n";
else output << "insert " << name() << ":RhoPWaveMagnitude " << ix << " "
<< _rhomagP[ix] << "\n";
}
for(unsigned int ix=0;ix<_rhophaseP.size();++ix) {
if(ix<2) output << "newdef " << name() << ":RhoPWavePhase " << ix << " "
<< _rhophaseP[ix] << "\n";
else output << "insert " << name() << ":RhoPWavePhase " << ix << " "
<< _rhophaseP[ix] << "\n";
}
for(unsigned int ix=0;ix<_rhomagD.size();++ix) {
if(ix<2) output << "newdef " << name() << ":RhoDWaveMagnitude " << ix << " "
<< _rhomagD[ix]*MeV2 << "\n";
else output << "insert " << name() << ":RhoDWaveMagnitude " << ix << " "
<< _rhomagD[ix]*MeV2 << "\n";
}
for(unsigned int ix=0;ix<_rhophaseD.size();++ix) {
if(ix<2) output << "newdef " << name() << ":RhoDWavePhase " << ix << " "
<< _rhophaseD[ix] << "\n";
else output << "insert " << name() << ":RhoDWavePhase " << ix << " "
<< _rhophaseD[ix] << "\n";
}
// use local values of the masses etc.
output << "newdef " << name() << ":LocalParameters " << _localparameters << "\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";
}
// maximum weights for the different channels
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;
}
diff --git a/Decay/VectorMeson/a1ThreePionCLEODecayer.h b/Decay/VectorMeson/a1ThreePionCLEODecayer.h
--- a/Decay/VectorMeson/a1ThreePionCLEODecayer.h
+++ b/Decay/VectorMeson/a1ThreePionCLEODecayer.h
@@ -1,604 +1,613 @@
// -*- C++ -*-
//
// a1ThreePionCLEODecayer.h 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.
//
#ifndef HERWIG_a1ThreePionCLEODecayer_H
#define HERWIG_a1ThreePionCLEODecayer_H
//
// This is the declaration of the a1ThreePionCLEODecayer class.
//
-#include "Herwig/Decay/DecayIntegrator.h"
-#include "Herwig/Decay/DecayPhaseSpaceMode.h"
+#include "Herwig/Decay/DecayIntegrator2.h"
+#include "Herwig/Decay/PhaseSpaceMode.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Helicity/LorentzPolarizationVector.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
*
* The <code>a1ThreePionCLEODecayer</code> class is designed to implement the decay
* of the \f$a_1\f$ to three pions using the model of Phys.Rev.D61:012002,2000,
* (hep-ex/9902022) (CLEO) which was fitted to the one charged and two neutral pion
* channel for the charged \f$a_1\f$ decay in \f$\tau \to a_1 -> \pi\pi\pi\f$.
* The other modes are infered from this using isospin. This is a sophisticated model
* including the coupling of the \f$a_1\f$ to the \f$\rho\f$, \f$\rho(1450)\f$,
* \f$f(1370)\f$ and \f$\sigma\f$ sigma mesons.
*
* In this case the current is given by
* \f[\mathcal{M} = \epsilon_\mu
* \left[F_1(p_2-p_3)^\mu+F_2(p_3-p_1)^\mu+F_3(p_1-p_2)^\mu\right].\f]
*
*
* The form factors for the \f$a_1^0 \to \pi^0 \pi^0 \pi^0\f$ mode are
* \f[F_1=
* \phantom{-}\frac23\left(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3)\right)
* -\frac23\left(g_\sigma B^S_\sigma(s_2)+g_{f_0}B^S_{f_0}(s_2)\right)
* +g_{f_2}\left(\frac12(s_3-s_2)B^D_{f_2}(s_1)
* -\frac1{18}\frac{(4m_{\pi^0}^2-s_2)(q^2+s_2-m_{\pi^0}^2)}{s_2}B^D_{f_2}(s_2)
* +\frac1{18}\frac{(4m_{\pi^0}^2-s_3)(q^2-m_{\pi^0}^2+s_3)}{s_3}B^D_{f_2}(s_3)\right)
*\f]
*
* \f[F_2=\phantom{-}\frac23(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3))
* -\frac23(g_\sigma B^S_\sigma(s_1)+g_{f_0}B^S_{f_0}(s_1))
* +g_{f_2}\left( \frac12(s_3-s1)B^D_{f_2}(s_2)
* -\frac1{18}\frac{(4m_{\pi^0}^2-s_1)(q^2+s_1-m_{\pi^0}^2)}{s_1}B^D_{f_2}(s_1)
* +\frac1{18}\frac{(4m_{\pi^0}^2-s_3)(q^2-m_{\pi^0}^2+s_3)}{s_3}B^D_{f_2}(s_3)\right)
*\f]
* \f[F_3=-\frac23(g_\sigma B^S_\sigma(s_1)+g_{f_0}B^S_{f_0}(s_1))
* +\frac23(g_\sigma B^S_\sigma(s_2)+g_{f_0}B^S_{f_0}(s_2))
* +g_{f_2}\left( \frac12(s_1-s_2)B^D_{f_2}(s_3)
* -\frac1{18}\frac{(4m_{\pi^0}^2-s_1)(q^2+s_1-m_{\pi^0}^2)}{s_1}B^D_{f_2}(s_1)
* +\frac1{18}\frac{(4m_{\pi^0}^2-s_2)(q^2+s_2-m_{\pi^0}^2)}{s_2}B^D_{f_2}(s_2)\right)
*\f]
*
* The form factors for the \f$a_1^+ \to \pi^0 \pi^0 \pi^+\f$ mode are
*
* \f[F_1=\sum_k\left\{-\frac{g^P_{\rho_k}}3B_{\rho_k}^P(s_1)
* -g^D_{\rho_k}B_{\rho_k}^P(s_2)
* \left((s_3-m_{\pi^+}^2)-(s_1-m_{\pi^0}^2)\right)\right\}
* +\frac23\left(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3)\right)
* +\frac{g_{f_2}}{18s_3}(q^2-m_{\pi^+}^2+s_3)(4m_{\pi^0}^2-s_3)B^D_{f_2}(s_3)
*\f]
*
* \f[F_2=\sum_k\left\{-\frac13g^P_{\rho_k}B_{\rho_k}^P(s_2)
* -g^D_{\rho_k}B_{\rho_k}^P(s_1)
* \left((s_3-m_{\pi^+}^2)-(s_2-m_{\pi^0}^2)\right)\right\}
* +\frac23\left(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3)\right)
* +\frac1{18s_3}g_{f_2}(q^2-m_{\pi^+}^2+s_3)(4m_{\pi^0}^2-s_3)B^D_{f_2}(s_3)
*\f]
*
* \f[F_3=\sum_k g^D_{\rho_k}\left\{
* -\frac13B_{\rho_k}^P(s_1)\left((s_3-m_{\pi^+}^2)-(s_2-m_{\pi^0}^2)\right)
* +\frac13B_{\rho_k}^P(s_2)\left((s_3-m_{\pi^+}^2)-(s_1-m_{\pi^0}^2)\right)\right\}
* -\frac{g_{f_2}}2(s_1-s_2)B^D_{f_2}(s_3)\f]
*
* The form factors for \f$a_1^0\to\pi^+\pi^-\pi^0\f$.
*
* \f[F_1=\sum_k\left\{g^P_{\rho_k}B_{\rho_k}^P(s_1)
* -\frac{g^D_{\rho_k}}3B_{\rho_k}^P(s_2)(s_3-m_{\pi^0}^2-s_1+m_{\pi^+}^2)\right\}
* +\frac23\left(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3)\right)
* +\frac{g_{f_2}}{18s_3}(q^2-m_{\pi^0}^2+s_3)(4m_{\pi^+}^2-s_3)B^D_{f_2}(s_3)\f]
*
* \f[F_2=\sum_k\left\{g^P_{\rho_k}B_{\rho_k}^P(s_2)
* -\frac{g^D_{\rho_k}}3B_{\rho_k}^P(s_1)(s_3-m_{\pi^0}^2-s_2+m_{\pi^+}^2)\right\}
* +\frac23\left(g_\sigma B^S_\sigma(s_3)+g_{f_0}B^S_{f_0}(s_3)\right)
* +\frac{g_{f_2}}{18s_3}(q^2-m_{\pi^0}^2+s_3)(4m_{\pi^+}^2-s_3)B^D_{f_2}(s_3)\f]
*
* \f[F_3=\sum_k
* g^D_{\rho_k}\left\{-\frac13B_{\rho_k}^P(s_1)(s_3-m_{\pi^0}^2-s_2+m_{\pi^+}^2)
* +\frac13B_{\rho_k}^P(s_2)(s_3-m_{\pi^0}^2-s_1+m_{\pi^+}^2)
* \right\}
* -\frac{g_{f_2}}2(s_1-s_2)B^D_{f_2}(s_3)\f]
*
* The form factors for \f$a_1^+\to \pi^+ \pi^+ \pi^-\f$ mode
*
* \f[F_1=\sum_k\left\{-g^P_{\rho_k}B_{\rho_k}^P(s_1)
* -\frac{g^D_{\rho_k}}3B_{\rho_k}^P(s_2)(s_1-s_3)\right\}
* -\frac23\left(g_\sigma B^S_\sigma(s_2)+g_{f_0} B^S_{f_0}(s_2)\right)
* +g_{f_2}\left(\frac12(s_3-s_2)B^D_{f_2}(s_1)
* -\frac1{18s_2}(4m_{\pi^+}^2-s_2)(q^2+s_2-m_{\pi^+}^2)B^D_{f_2}(s_2)\right)\f]
*
* \f[F_2=\sum_k\left\{-g^P_{\rho_k}B_{\rho_k}^P(s_2)
* -\frac{g^D_{\rho_k}}3B_{\rho_k}^P(s_1)(s_2-s_3)\right\}
* -\frac23\left(g_\sigma B^S_\sigma(s_1)+g_{f_0} B^S_{f_0}(s_1)\right)
* +g_{f_2}\left(\frac12(s_3-s_1)B^D_{f_2}(s_2)
* -\frac1{18s_1}(4m_{\pi^+}^2-s_1)(q^2+s_1-m_{\pi^+}^2)B^D_{f_2}(s_1)\right)\f]
*
* \f[F_3=\sum_k
* -g^D_{\rho_k}\left( \frac13(s_2-s_3)B_{\rho_k}^P(s_1)
* -\frac13(s_1-s_3)B_{\rho_k}^P(s_2)\right)
* -\frac23\left(g_\sigma B^S_\sigma(s_1)+g_{f_0}B^S_{f_0}(s_1)\right)
* +\frac23\left(g_\sigma B^S_\sigma(s_2)+g_{f_0}B^S_{f_0}(s_2)\right)\f]
*\f[
* +g_{f_2}\left(-\frac1{18s_1}(4m_{\pi^+}^2-s_1)(q^2+s_1-m_{\pi^+}^2)B^D_{f_2}(s_1)
* +\frac1{18s_2}(4m_{\pi^+}^2-s_2)(q^2+s_2-m_{\pi^+}^2)B^D_{f_2}(s_2)\right)\f]
*
* where
*
* - \f$g_{f_2}\f$ is the coupling of the \f$f_2\f$ to the \f$a_1\f$
* - \f$g_{f_0}\f$ is the coupling of the \f$f_0(1370)\f$ to the \f$a_1\f$
* - \f$g_{\sigma}\f$ is the coupling of the \f$\sigma\f$ to the \f$a_1\f$
* - \f$g^P_{\rho_k}\f$ is the \f$p\f$-wave coupling of the \f$\rho_k\f$ multiplet
* to the \f$a_1\f$.
* - \f$g^D_{\rho_k}\f$ is the \f$d\f$-wave coupling of the \f$\rho_k\f$ multiplet
* to the \f$a_1\f$.
* - \f$s_3=m^2_{12}\f$ is the invariant mass squared of particles 1 and 2.
* - \f$s_2=m^2_{13}\f$ is the invariant mass squared of particles 1 and 3.
* - \f$s_1=m^2_{23}\f$ is the invariant mass squared of particles 2 and 3.
*
* The Breit-Wigner factors are given by
\f$B^L_Y(s_i) = \frac{m^2_Y}{m^2_Y-s_i)+im_Y\Gamma^{Y,L}(s_i)}\f$
* where
* \f$\Gamma^{Y,L}(s_i) = \Gamma^Y\left(\frac{p(s_i)}{p(M_Y}\right)^{2L+1}\frac{m_Y}{\sqrt{s_i}}\f$
* \f$m_Y\f$ and \f$\Gamma^Y\f$ are the mass and width of the particle \f$Y\f$
* respectively. \f$p(s_i)\f$ is the momentum of the outgoing pion in the
* rest frame of the resonanc \f$Y\f$.
*
* @see ThreePionCLEOCurrent
- * @see DecayIntegrator
+ * @see DecayIntegrator2
*
*/
-class a1ThreePionCLEODecayer: public DecayIntegrator {
+class a1ThreePionCLEODecayer: public DecayIntegrator2 {
public:
/**
* Default constructor.
*/
a1ThreePionCLEODecayer();
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const;
-
+
/**
* 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 decay The particles produced in the decay.
+ * @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 ParticleVector & decay, MEOption meopt) const;
+ 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;
/**
* Method to return an object to calculate the 3 body partial width.
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
/**
* The matrix element to be integrated for the three-body decays as a function
* of the invariant masses of pairs of the outgoing particles.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The matrix element
*/
virtual double threeBodyMatrixElement(const int imode , const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy m1,
const Energy m2, const Energy m3) const;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) 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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const { return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const { return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving and
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object to the begining of the run phase.
*/
virtual void doinitrun();
//@}
private:
/**
* Private and non-existent assignment operator.
*/
a1ThreePionCLEODecayer & operator=(const a1ThreePionCLEODecayer &);
private:
/**
* Breit wigner for the \f$\rho\f$, \f$B^P_{\rho_k}(q^2)\f$.
* @param ires The \f$\rho\f$ multiplet to used.
* @param q2 The scale, \f$q^2\f$.
* @param icharge Which pion masses to use for the momentum calculation
* @return The Breit-Wigner
*/
Complex rhoBreitWigner(int ires, Energy2 q2,int icharge) const {
Energy q=sqrt(q2);
Complex ii(0.,1.);
double ratio = icharge==0 ?
Kinematics::pstarTwoBodyDecay(q,_mpic,_mpic)/_prhocc[ires] :
Kinematics::pstarTwoBodyDecay(q,_mpic,_mpi0)/_prhoc0[ires];
Energy gamrun=_rhowidth[ires]*pow(ratio,3)*_rhomass[ires]/q;
return sqr(_rhomass[ires])
/(sqr(_rhomass[ires])-q2-ii*_rhomass[ires]*gamrun);
}
/**
* Breit wigner for the \f$\sigma\f$, \f$B^S_\sigma(q^2)\f$.
* @param q2 The scale, \f$q^2\f$.
* @param icharge Which pion masses to use for the momentum calculation
* @return The Breit-Wigner
*/
Complex sigmaBreitWigner(Energy2 q2,int icharge) const {
Energy q=sqrt(q2);
Complex ii(0.,1.);
double ratio = icharge==0 ?
Kinematics::pstarTwoBodyDecay(q,_mpic,_mpic)/_psigmacc :
Kinematics::pstarTwoBodyDecay(q,_mpi0,_mpi0)/_psigma00;
Energy gamrun=_sigmawidth*ratio*_sigmamass/q;
return sqr(_sigmamass)/(sqr(_sigmamass)-q2-ii*_sigmamass*gamrun);
}
/**
* Breit wigner for the \f$f_0(1370)\f$, \f$B^S_{f_0}(q^2)\f$.
* @param q2 The scale, \f$q^2\f$.
* @param icharge Which pion masses to use for the momentum calculation
* @return The Breit-Wigner
*/
Complex f0BreitWigner(Energy2 q2,int icharge) const {
Energy q=sqrt(q2);
Complex ii(0.,1.);
double ratio = icharge==0 ?
Kinematics::pstarTwoBodyDecay(q,_mpic,_mpic)/_pf0cc :
Kinematics::pstarTwoBodyDecay(q,_mpi0,_mpi0)/_pf000;
Energy gamrun=_f0width*ratio*_f0mass/q;
return sqr(_f0mass)/(sqr(_f0mass)-q2-ii*_f0mass*gamrun);
}
/**
* Breit wigner for the \f$f_2\f$, \f$B^D_{f_2}(q^2)\f$.
* @param q2 The scale, \f$q^2\f$.
* @param icharge Which pion masses to use for the momentum calculation
* @return The Breit-Wigner
*/
Complex f2BreitWigner(Energy2 q2,int icharge) const {
Energy q=sqrt(q2);
Complex ii(0.,1.);
double ratio = icharge==0 ?
Kinematics::pstarTwoBodyDecay(q,_mpic,_mpic)/_pf2cc :
Kinematics::pstarTwoBodyDecay(q,_mpi0,_mpi0)/_pf200;
Energy gamrun=_f2width*pow(ratio,5)*_f2mass/q;
return sqr(_f2mass)/(sqr(_f2mass)-q2-ii*_f2mass*gamrun);
}
/**
* Calculate the form factors
* @param iopt The mode being calculated in the order given above
* @param ichan The phase space channel in the order given in the doinit member.
* @param q2 The sacale \f$q^2\f$.
* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
* @param F1 The form factor \f$F_1\f$.
* @param F2 The form factor \f$F_2\f$.
* @param F3 The form factor \f$F_3\f$.
*
*/
void formFactors(int iopt,int ichan,Energy2 q2,Energy2 s1,Energy2 s2,
Energy2 s3,
complex<InvEnergy> & F1,
complex<InvEnergy> & F2,
complex<InvEnergy> & F3) const;
private:
/**
* Masses of the rho resonaces
*/
vector<Energy> _rhomass;
/**
* Widths of the rho resonaces
*/
vector<Energy> _rhowidth;
/**
* Momentum of the particles produced in charged rho decay
*/
vector<Energy> _prhocc;
/**
* Momentum of the particles produced in neutral rho decay
*/
vector<Energy> _prhoc0;
/**
* Mass of the \f$f_2\f$.
*/
Energy _f2mass;
/**
* Width of the \f$f_2\f$.
*/
Energy _f2width;
/**
* Momentum for the decay of the \f$f_2\f$ to two charged pions.
*/
Energy _pf2cc;
/**
* Momentum for the decay of the \f$f_2\f$ to two neutral pions.
*/
Energy _pf200;
/**
* Mass of the \f$f_0(1370)\f$.
*/
Energy _f0mass;
/**
* Width of the \f$f_0(1370)\f$.
*/
Energy _f0width;
/**
* Momentum for the decay of the \f$f_0(1370)\f$ to two charged pions.
*/
Energy _pf0cc;
/**
* Momentum for the decay of the \f$f_0(1370)\f$ to two neutral pions.
*/
Energy _pf000;
/**
* Mass of the \f$\sigma\f$ meson.
*/
Energy _sigmamass;
/**
* Width of the \f$\sigma\f$ meson.
*/
Energy _sigmawidth;
/**
* Momentum for the decay of the \f$\sigma\f$ to two charged pions.
*/
Energy _psigmacc;
/**
* Momentum for the decay of the \f$\sigma\f$ to two neutral pions.
*/
Energy _psigma00;
/**
* Mass of the neutral pion
*/
Energy _mpi0;
/**
* Mass of the charged pion
*/
Energy _mpic;
/**
* overall coupling for the decay
*/
InvEnergy _coupling;
/**
* Magnitude of the \f$p\f$-wave couplings of the rho resonance, \f$g^P_{\rho_k}\f$,
* (\f$\beta_{1,2}\f$ in the CLEO paper.)
*/
vector<double> _rhomagP;
/**
* Phase of the \f$p\f$-wave couplings of the rho resonance, \f$g^P_{\rho_k}\f$,
* (\f$\beta_{1,2}\f$ in the CLEO paper.)
*/
vector<double> _rhophaseP;
/**
*\f$p\f$-wave couplings of the rho resonance, \f$g^P_{\rho_k}\f$,
* (\f$\beta_{1,2}\f$ in the CLEO paper.)
*/
vector<Complex> _rhocoupP;
/**
* Magnitude of the \f$d\f$-wave couplings of the rho resonance, \f$g^D_{\rho_k}\f$,
* (\f$\beta_{3,4}\f$ in the CLEO paper.)
*/
vector<InvEnergy2> _rhomagD;
/**
* Phase of the \f$d\f$-wave couplings of the rho resonance, \f$g^D_{\rho_k}\f$,
* (\f$\beta_{3,4}\f$ in the CLEO paper.)
*/
vector<double>_rhophaseD;
/**
* \f$d\f$-wave couplings of the rho resonance, \f$g^D_{\rho_k}\f$,
* (\f$\beta_{3,4}\f$ in the CLEO paper.)
*/
vector<complex<InvEnergy2> > _rhocoupD;
/**
* Magntiude of the coupling of the \f$f_2\f$ resonance, \f$g_{f_2}\f$,
* (\f$\beta_5\f$ in the CLEO paper.)
*/
InvEnergy2 _f2mag;
/**
* Phase of the coupling of the \f$f_2\f$ resonance, \f$g_{f_2}\f$,
* (\f$\beta_5\f$ in the CLEO paper.)
*/
double _f2phase;
/**
* Coupling of the \f$f_2\f$ resonance, \f$g_{f_2}\f$,
* (\f$\beta_5\f$ in the CLEO paper.)
*/
complex<InvEnergy2> _f2coup;
/**
* Magntiude of the coupling of the \f$f_0(1370)\f$ resonance, \f$g_{f_0}\f$,
* (\f$\beta_6\f$ in the CLEO paper.)
*/
double _f0mag;
/**
* Phase of the coupling of the \f$f_0(1370)\f$ resonance, \f$g_{f_0}\f$,
* (\f$\beta_6\f$ in the CLEO paper.)
*/
double _f0phase;
/**
* Coupling of the \f$f_0(1370)\f$ resonance, \f$g_{f_0}\f$,
* (\f$\beta_6\f$ in the CLEO paper.)
*/
Complex _f0coup;
/**
* Magntiude of the coupling of the \f$\sigma\f$ resonance, \f$g_\sigma\f$,
* (\f$\beta_7\f$ in the CLEO paper.)
*/
double _sigmamag;
/**
* Phase of the coupling of the \f$\sigma\f$ resonance, \f$g_\sigma\f$,
* (\f$\beta_7\f$ in the CLEO paper.)
*/
double _sigmaphase;
/**
* Coupling of the \f$\sigma\f$ resonance, \f$g_\sigma\f$,
* (\f$\beta_7\f$ in the CLEO paper.)
*/
Complex _sigmacoup;
/**
* Use local values of the mass parameters
*/
bool _localparameters;
/**
* Weights for the channels for the zero charged pion channel.
*/
mutable vector<double> _zerowgts;
/**
* Weights for the channels for the one charged pion channel.
*/
mutable vector<double> _onewgts;
/**
* Weights for the channels for the two charged pion channel.
*/
mutable vector<double> _twowgts;
/**
* Weights for the channels for the three charged pion channel.
*/
mutable vector<double> _threewgts;
/**
* Maximum weight for the zero charged pion channel.
*/
mutable double _zeromax;
/**
* Maximum weight for the one charged pion channel.
*/
mutable double _onemax;
/**
* Maximum weight for the two charged pion channel.
*/
mutable double _twomax;
/**
* Maximum weight for the three charged pion channel.
*/
mutable double _threemax;
/**
* Spin density matrix
*/
mutable RhoDMatrix _rho;
/**
* Polarization vectors
*/
mutable vector<Helicity::LorentzPolarizationVector> _vectors;
};
}
#endif /* HERWIG_a1ThreePionCLEODecayer_H */
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Mon, Jan 20, 9:26 PM (1 d, 1 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242470
Default Alt Text
(69 KB)
Attached To
rHERWIGHG herwighg
Event Timeline
Log In to Comment