Page MenuHomeHEPForge
Files F7879518

No OneTemporary
Actions

View Options

diff --git a/Decay/General/VectorCurrentDecayer.cc b/Decay/General/VectorCurrentDecayer.cc
--- a/Decay/General/VectorCurrentDecayer.cc
+++ b/Decay/General/VectorCurrentDecayer.cc
@@ -1,273 +1,273 @@
// -*- C++ -*-
//
// VectorCurrentDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 VectorCurrentDecayer class.
//
#include "VectorCurrentDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "Herwig/Models/General/BSMModel.h"
using namespace Herwig;
IBPtr VectorCurrentDecayer::clone() const {
return new_ptr(*this);
}
IBPtr VectorCurrentDecayer::fullclone() const {
return new_ptr(*this);
}
void VectorCurrentDecayer::persistentOutput(PersistentOStream & os) const {
os << inpart_ << currentOut_ << current_ << mode_ << wgtloc_ << wgtmax_ << weights_ << cSMmed_;
}
void VectorCurrentDecayer::persistentInput(PersistentIStream & is, int) {
is >> inpart_ >> currentOut_ >> current_ >> mode_ >> wgtloc_ >> wgtmax_ >> weights_ >> cSMmed_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<VectorCurrentDecayer,DecayIntegrator>
describeHerwigVectorCurrentDecayer("Herwig::VectorCurrentDecayer", "Herwig.so");
void VectorCurrentDecayer::Init() {
static ClassDocumentation<VectorCurrentDecayer> documentation
("The VectorCurrentDecayer class is designed for the decays of low mass vector bosons");
}
void VectorCurrentDecayer::setDecayInfo(PDPtr in, const vector<tPDPtr> & outCurrent, WeakCurrentPtr current) {
inpart_ = in;
currentOut_ = outCurrent;
current_ = current;
// cast the model
Ptr<BSMModel>::ptr model = dynamic_ptr_cast<Ptr<BSMModel>::ptr>(generator()->standardModel());
bool foundU(false),foundD(false),foundS(false);
// find the vertices we need and extract the couplings
for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) {
VertexBasePtr vertex = model->vertex(ix);
if(vertex->getNpoint()!=3) continue;
for(unsigned int iloc = 0;iloc < 3; ++iloc) {
vector<long> ext = vertex->search(iloc, in->id());
if(ext.empty()) continue;
for(unsigned int ioff=0;ioff<ext.size();ioff+=3) {
if(iloc!=2) assert(false);
if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) {
foundD = true;
vertex->setCoupling(sqr(in->mass()),getParticleData(1),getParticleData(-1),in);
cSMmed_[0] = vertex->norm();
}
else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) {
foundU = true;
vertex->setCoupling(sqr(in->mass()),getParticleData(2),getParticleData(-2),in);
cSMmed_[1] = vertex->norm();
}
else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) {
foundS = true;
vertex->setCoupling(sqr(in->mass()),getParticleData(3),getParticleData(-3),in);
cSMmed_[2] = vertex->norm();
}
}
}
}
if(!foundD) {
throw InitException() << "Cannot find down quark coupling in VectorCurrentDecayer::doinit()";
}
if(!foundU) {
throw InitException() << "Cannot find up quark coupling in VectorCurrentDecayer::doinit()";
}
if(!foundS) {
throw InitException() << "Cannot find strange quark coupling in VectorCurrentDecayer::doinit()";
}
}
int VectorCurrentDecayer::modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const {
vector<long> id;
id.push_back(parent->id());
for(unsigned int ix=0;ix<children.size();++ix) id.push_back(children[ix]->id());
return modeNumber(cc,id);
}
void VectorCurrentDecayer::doinitrun() {
current_->initrun();
DecayIntegrator::doinitrun();
}
void VectorCurrentDecayer::doinit() {
DecayIntegrator::doinit();
// make sure the current got initialised
current_->init();
// find the mode
for(unsigned int ix=0;ix<current_->numberOfModes();++ix) {
// get the external particles for this mode
int iq(0),ia(0);
tPDVector ptemp = current_->particles(inpart_->iCharge(),ix,iq,ia);
// check this is the right mode
if(ptemp.size()!=currentOut_.size()) continue;
vector<bool> matched(ptemp.size(),false);
bool match = true;
for(unsigned int iy=0;iy<currentOut_.size();++iy) {
bool found = false;
for(unsigned int iz=0;iz<ptemp.size();++iz) {
if(!matched[iz]&&ptemp[iz]==currentOut_[iy]) {
found = true;
matched[iz] = true;
break;
}
}
if(!found) {
match = false;
break;
}
}
if(!match) continue;
tPDVector out = {};
out.insert(std::end(out), std::begin(ptemp), std::end(ptemp));
// create the mode
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(inpart_,out,1.));
PhaseSpaceChannel channel(mode,true);
bool done=current_->createMode(inpart_->iCharge(),tcPDPtr(),FlavourInfo(),
ix,mode,0,-1,channel,inpart_->mass());
if(done) {
// the maximum weight and the channel weights
// the weights for the channel
if(weights_.empty()) {
weights_.resize(mode->channels().size(),1./(mode->channels().size()));
}
mode_ = ix;
// special for the two body modes
if(out.size()==2) {
weights_.clear();
mode=new_ptr(PhaseSpaceMode(inpart_,out,1.));
}
mode->maxWeight(wgtmax_);
mode->setWeights(weights_);
addMode(mode);
}
break;
}
}
int VectorCurrentDecayer::modeNumber(bool & cc, vector<long> id) const {
// incoming particle
long idtemp;
tPDPtr p0=getParticleData(id[0]);
idtemp = p0->CC() ? -id[0] : id[0];
if( id[0] ==inpart_->id()) cc=false;
else if(idtemp==inpart_->id()) cc=true ;
else return -1;
vector<int> idout;
for(vector<long>::iterator it=++id.begin();it!=id.end();++it) {
idout.push_back(*it);
}
unsigned int icurr=current_->decayMode(idout);
if(mode_==icurr) return 0;
else return -1;
}
void VectorCurrentDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming,true,false);
weakCurrent()->constructSpinInfo(ParticleVector(decay.begin(),decay.end()));
}
double VectorCurrentDecayer::me2(const int ichan, const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
using namespace ThePEG::Helicity;
// polarization vectors for the incoming particle
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
const_ptr_cast<tPPtr>(&part),
incoming,false);
// fix rho if no correlations
fixRho(rho_);
}
// work out the mapping for the hadron vector
int nOut = momenta.size();
vector<unsigned int> constants(nOut+1);
vector<PDT::Spin > iSpin(nOut);
vector<int> hadrons(nOut);
int itemp(1);
int ix(nOut);
do {
--ix;
iSpin[ix] = outgoing[ix]->iSpin();
itemp *= iSpin[ix];
constants[ix] = itemp;
hadrons[ix] = outgoing[ix]->id();
}
while(ix>0);
constants[nOut] = 1;
Energy2 scale(sqr(part.mass()));
// calculate the hadron current
Energy q = part.mass();
// currents for the different flavour components
vector<LorentzPolarizationVectorE>
hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero),
mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero),
mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar),
mode(),ichan,q,outgoing,momenta,DecayIntegrator::Calculate));
// compute the matrix element
GeneralDecayMEPtr newME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,iSpin)));
vector<unsigned int> ihel(momenta.size()+1);
unsigned int hI0_size = hadronI0.size();
unsigned int hI1_size = hadronI1.size();
unsigned int hss_size = hadronssbar.size();
unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size);
for(unsigned int hhel=0;hhel<maxsize;++hhel) {
// map the index for the hadrons to a helicity state
for(int ix=nOut;ix>0;--ix) {
ihel[ix]=(hhel%constants[ix-1])/constants[ix];
}
for(ihel[0]=0;ihel[0]<3;++ihel[0]) {
Complex amp = 0.;
// work on coefficients for the I1 and I0 bits
if(hI0_size != 0 )
- amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel]));
+ amp += Complex((cSMmed_[0]+cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI0[hhel])));
if(hI1_size !=0)
- amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel]));
+ amp += Complex((cSMmed_[0]-cSMmed_[1])/sqrt(2.)/q*(vectors_[ihel[0]].wave().dot(hadronI1[hhel])));
if(hss_size !=0)
- amp += cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel]));
+ amp += Complex(cSMmed_[2] /q*(vectors_[ihel[0]].wave().dot(hadronssbar[hhel])));
(*newME)(ihel) = amp;
}
}
// store the matrix element
ME(newME);
// return the answer
double output = (ME()->contract(rho_)).real();
return output;
}
Energy VectorCurrentDecayer::partialWidth(tPDPtr part, vector<tPDPtr> out) {
vector<long> id;
id.push_back(part->id());
for(unsigned int ix=0;ix<out.size();++ix) id.push_back(out[ix]->id());
bool cc;
int mode=modeNumber(cc,id);
imode(mode);
return initializePhaseSpaceMode(mode,true,true);
}
diff --git a/Decay/ScalarMeson/PScalar4FermionsDecayer.cc b/Decay/ScalarMeson/PScalar4FermionsDecayer.cc
--- a/Decay/ScalarMeson/PScalar4FermionsDecayer.cc
+++ b/Decay/ScalarMeson/PScalar4FermionsDecayer.cc
@@ -1,398 +1,398 @@
// -*- C++ -*-
//
// PScalar4FermionsDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 PScalar4FermionsDecayer class.
//
#include "PScalar4FermionsDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Helicity/ScalarSpinInfo.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void PScalar4FermionsDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<_incoming.size();++ix)
_maxweight[ix] = mode(ix)->maxWeight();
}
}
PScalar4FermionsDecayer::PScalar4FermionsDecayer()
: _coupling(1,0.025159062/GeV), _incoming(1,111), _outgoing1(1,11),
_outgoing2(1,11), _maxweight(1,0.000234211),
_includeVMD(1,2),_VMDid(1,113), _VMDmass(1,0.7758*GeV),
_VMDwidth(1,0.1503*GeV), _initsize(1) {
// intermediates
generateIntermediates(false);
}
void PScalar4FermionsDecayer::doinit() {
DecayIntegrator::doinit();
// check the parameters are consistent
unsigned int isize=_coupling.size();
if(isize!=_incoming.size() || isize!=_outgoing1.size() || isize!=_outgoing2.size()||
isize!=_maxweight.size() || isize!=_includeVMD.size()|| isize!=_VMDid.size() ||
isize!=_VMDmass.size() || isize!=_VMDwidth.size())
throw InitException() << "Inconsistent parameters in PScalar4FermionsDecayer"
<< Exception::abortnow;
// create the integration channels for each mode
tPDPtr gamma=getParticleData(ParticleID::gamma);
for(unsigned int ix=0;ix<_incoming.size();++ix) {
tPDPtr in = getParticleData(_incoming[ix]);
tPDVector out={getParticleData( _outgoing1[ix]),
getParticleData(-_outgoing1[ix]),
getParticleData( _outgoing2[ix]),
getParticleData(-_outgoing2[ix])};
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix]));
PhaseSpaceChannel newChannel((PhaseSpaceChannel(mode),0,gamma,0,gamma,1,1,1,2,2,3,2,4));
newChannel.setJacobian(1,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1);
newChannel.setJacobian(2,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1);
mode->addChannel(newChannel);
PhaseSpaceChannel newChannel2((PhaseSpaceChannel(mode),0,gamma,0,gamma,1,3,1,2,2,1,2,4));
newChannel2.setJacobian(1,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1);
newChannel2.setJacobian(2,PhaseSpaceChannel::PhaseSpaceResonance::Power,-1.1);
mode->addChannel(newChannel2);
addMode(mode);
}
// set up the values for the VMD factor if needed (copy the default mass and width
// into the array)
for(unsigned ix=0;ix<isize;++ix) {
if(_includeVMD[ix]==1) {
_VMDmass[ix]=getParticleData(_VMDid[ix])->mass();
_VMDwidth[ix]=getParticleData(_VMDid[ix])->width();
}
}
}
int PScalar4FermionsDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
// must be four outgoing particles
if(children.size()!=4) return -1;
// get the id's of the outgoing particles
int id[4]={0,0,0,0};
bool done[4]={false,false,false,false};
unsigned int ix(0),iy(0);
// ids of the particles
int id0(parent->id()),idtemp(-1),idl1(-1),idl2(-1),idt[2];
tPDVector::const_iterator pit = children.begin();
for ( ;pit!=children.end();++pit) {
id[ix]=(**pit).id();
done[ix]=false;
++ix;
}
// find the two lepton pairs
// find the first fermion
ix=0;
do {
if( id[ix]>0 && !done[ix] ) {
done[ix]=true;
idtemp=id[ix];
}
++ix;
}
while(ix<4&&idtemp<0);
if(idtemp<0) return -1;
// find its antiparticle
ix=0;
do {
if( id[ix]==-idtemp && !done[ix] ) {
done[ix]=true;
idl1=idtemp;
}
++ix;
} while( ix<4 && idl1<0 );
if(idl1<0) return -1;
// find the second particle antiparticle pair
for(ix=0;ix<4;++ix) {
if(!done[ix]) {
idt[iy]=id[ix];
++iy;
}
}
if(idt[0]==-idt[1]) idl2=abs(idt[0]);
if(idl2<0) return -1;
// loop over the modes and see if this is one of them
ix=0;
int imode(-1);
do {
if(_incoming[ix]==id0) {
if((idl1==_outgoing1[ix]&&idl2==_outgoing2[ix])||
(idl2==_outgoing1[ix]&&idl1==_outgoing2[ix])) imode=ix;
}
++ix;
}
while(imode<0&&ix<_incoming.size());
cc=false;
return imode;
}
void PScalar4FermionsDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_coupling,1/MeV)
<< _incoming << _outgoing1 << _outgoing2 << _maxweight
<< _includeVMD << _VMDid
<< ounit(_VMDmass,MeV) << ounit(_VMDwidth,MeV);
}
void PScalar4FermionsDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_coupling,1/MeV)
>> _incoming >> _outgoing1 >> _outgoing2 >> _maxweight
>> _includeVMD >> _VMDid
>> iunit(_VMDmass,MeV) >> iunit(_VMDwidth,MeV);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PScalar4FermionsDecayer,DecayIntegrator>
describeHerwigPScalar4FermionsDecayer("Herwig::PScalar4FermionsDecayer", "HwSMDecay.so");
void PScalar4FermionsDecayer::Init() {
static ClassDocumentation<PScalar4FermionsDecayer> documentation
("The PScalar4FermionsDecayer class is designed for the decay"
" of a pseudosclar meson to four fermions. It is intended for the decay of"
"the pion to two electron-positron pairs.");
static ParVector<PScalar4FermionsDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&PScalar4FermionsDecayer::_incoming,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<PScalar4FermionsDecayer,int> interfaceOutcoming1
("Outgoing1",
"The PDG code for the first outgoing fermion",
&PScalar4FermionsDecayer::_outgoing1,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<PScalar4FermionsDecayer,int> interfaceOutcoming2
("Outgoing2",
"The PDG code for the second outgoing fermion",
&PScalar4FermionsDecayer::_outgoing2,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<PScalar4FermionsDecayer,InvEnergy> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&PScalar4FermionsDecayer::_coupling,
1/MeV, 0, ZERO, -10000/MeV, 10000/MeV, false, false, true);
static ParVector<PScalar4FermionsDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&PScalar4FermionsDecayer::_maxweight,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<PScalar4FermionsDecayer,int> interfaceIncludeVMD
("IncludeVMD",
"There are three options for 0 the VMD factor is not included, for 1 the factor "
"is included using the default mass and width of the particle specified by"
" VMDID, and for 2 the factor is included using the mass and width specified"
" by VMDwidth and VMDmass.",
&PScalar4FermionsDecayer::_includeVMD,
0, 0, 0, 0, 2, false, false, true);
static ParVector<PScalar4FermionsDecayer,int> interfaceVMDID
("VMDID",
"The PDG code for the particle to be used for the VMD factor.",
&PScalar4FermionsDecayer::_VMDid,
0, 0, 0, -10000, 10000, false, false, true);
static ParVector<PScalar4FermionsDecayer,Energy> interfaceVMDmass
("VMDmass",
"The mass to use for the particle in the VMD factor",
&PScalar4FermionsDecayer::_VMDmass,
1.*MeV, -1, ZERO, -10000*MeV, 10000*MeV, false, false, true);
static ParVector<PScalar4FermionsDecayer,Energy> interfaceVMDwidth
("VMDwidth",
"The width to use for the particle in the VMD factor",
&PScalar4FermionsDecayer::_VMDwidth,
1.*MeV, -1, ZERO, -10000*MeV, 10000*MeV, false, false, true);
}
void PScalar4FermionsDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
// set up the spin information for the decay products
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
// set up the spin information for the decay products
for(unsigned int ix=0;ix<2;++ix) {
SpinorBarWaveFunction::
constructSpinInfo(_wavebar[ix],decay[2*ix ],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(_wave[ix] ,decay[2*ix+1],outgoing,true);
}
}
double PScalar4FermionsDecayer::me2(const int,const Particle & part,
const tPDVector &,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half)));
bool identical((_outgoing1[imode()]==_outgoing2[imode()]));
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
}
// calculate the spinors
for(unsigned int ix=0;ix<2;++ix) {
_wavebar[ix].resize(2);
_wave [ix].resize(2);
for(unsigned int ihel=0;ihel<2;++ihel) {
_wavebar[ix][ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[2*ix ],ihel,Helicity::outgoing);
_wave [ix][ihel] = HelicityFunctions::dimensionedSpinor (-momenta[2*ix+1],ihel,Helicity::outgoing);
}
}
// // momenta of the outgoing photons
Lorentz5Momentum poff[4];
poff[0]=momenta[0]+momenta[1];
poff[0].rescaleMass();
poff[1]=momenta[2]+momenta[3];
poff[1].rescaleMass();
if(identical) {
poff[2]=momenta[2]+momenta[1];
poff[2].rescaleMass();
poff[3]=momenta[0]+momenta[3];
poff[3].rescaleMass();
}
// compute the currents for the two leptonic decays
LorentzPolarizationVectorE current[4][2][2];
unsigned int it,ix,iy,iz;
for(iz=0;iz<2;++iz) {
it = iz==0 ? 1 : 0;
for(ix=0;ix<2;++ix) {
for(iy=0;iy<2;++iy) {
current[iz][iy][ix] = _wave[iz][ix].vectorCurrent(_wavebar[iz][iy]);
// the second two currents
if(identical)
current[iz+2][iy][ix] = _wave[it][ix].vectorCurrent(_wavebar[iz][iy]);
}
}
}
// invariants
Energy2 m12(sqr(poff[0].mass()));
Energy2 m34(sqr(poff[1].mass()));
Energy2 m14(ZERO), m23(ZERO);
complex<InvEnergy4> prop1(1./m12/m34),prop2(0./sqr(MeV2));
Complex ii(0.,1.);
if(identical) {
m14=poff[2].mass()*poff[2].mass();
m23=poff[3].mass()*poff[3].mass();
prop2=1./m14/m23;
}
// the VMD factor if needed
if(_includeVMD[imode()]>0) {
Energy2 mrho2(_VMDmass[imode()]*_VMDmass[imode()]);
Energy2 mwrho(_VMDmass[imode()]*_VMDwidth[imode()]);
prop1 = prop1*(-mrho2+ii*mwrho)/(m12-mrho2+ii*mwrho)*
(-mrho2+ii*mwrho)/(m34-mrho2+ii*mwrho);
if(identical) {
prop2 = prop2*(-mrho2+ii*mwrho)/(m14-mrho2+ii*mwrho)*
(-mrho2+ii*mwrho)/(m23-mrho2+ii*mwrho);
}
}
// prefactor
Complex pre(_coupling[imode()]*4.*Constants::pi
*SM().alphaEM()*part.mass());
Complex diag;
// now compute the matrix element
LorentzVector<complex<Energy3> > eps;
vector<unsigned int> ispin(5,0);
for(ispin[1]=0; ispin[1]<2;++ispin[1]) {
for(ispin[2]=0;ispin[2]<2;++ispin[2]) {
for(ispin[3]=0;ispin[3]<2;++ispin[3]) {
for(ispin[4]=0;ispin[4]<2;++ispin[4]) {
// the first diagram
eps = epsilon(current[0][ispin[1]][ispin[2]],poff[1],
current[1][ispin[3]][ispin[4]]);
- diag = prop1*(eps*poff[0]);
+ diag = Complex(prop1*(eps*poff[0]));
// exchanged diagram if identical particles
// (sign due normal ordering)
if(identical) {
eps = epsilon(current[2][ispin[1]][ispin[4]],poff[3],
current[3][ispin[3]][ispin[2]]);
- diag-= prop2*(eps*poff[2]);
+ diag-= Complex(prop2*(eps*poff[2]));
}
(*ME())(ispin)=pre*diag;
}
}
}
}
double me=ME()->contract(_rho).real();
if(identical) me *= 0.25;
return me;
}
// output the setup info for the particle database
void PScalar4FermionsDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
// parameters for the DecayIntegrator base class
DecayIntegrator::dataBaseOutput(output,false);
for(unsigned int ix=0;ix<_incoming.size();++ix) {
if(ix<_initsize) {
output << "newdef " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "newdef " << name() << ":Outgoing1 " << ix << " "
<< _outgoing1[ix] << "\n";
output << "newdef " << name() << ":Outgoing2 " << ix << " "
<< _outgoing2[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*MeV << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
output << "newdef " << name() << ":IncludeVMD " << ix << " "
<< _includeVMD[ix] << "\n";
output << "newdef " << name() << ":VMDID " << ix << " "
<< _VMDid[ix] << "\n";
output << "newdef " << name() << ":VMDmass " << ix << " "
<< _VMDmass[ix]/MeV << "\n";
output << "newdef " << name() << ":VMDwidth " << ix << " "
<< _VMDwidth[ix]/MeV << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "insert " << name() << ":Outgoing1 " << ix << " "
<< _outgoing1[ix] << "\n";
output << "insert " << name() << ":Outgoing2 " << ix << " "
<< _outgoing2[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*MeV << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
output << "insert " << name() << ":IncludeVMD " << ix << " "
<< _includeVMD[ix] << "\n";
output << "insert " << name() << ":VMDID " << ix << " "
<< _VMDid[ix] << "\n";
output << "insert " << name() << ":VMDmass " << ix << " "
<< _VMDmass[ix]/MeV << "\n";
output << "insert " << name() << ":VMDwidth " << ix << " "
<< _VMDwidth[ix]/MeV << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";" << endl;
}
diff --git a/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc b/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc
--- a/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc
+++ b/Decay/ScalarMeson/PScalarPScalarVectorDecayer.cc
@@ -1,283 +1,283 @@
// -*- C++ -*-
//
// PScalarPScalarVectorDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 PScalarPScalarVectorDecayer class.
//
#include "PScalarPScalarVectorDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void PScalarPScalarVectorDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<_incoming.size();++ix)
if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight();
}
}
PScalarPScalarVectorDecayer::PScalarPScalarVectorDecayer() {
_incoming = { 100111, 100211, 100211, 100311, 100321, 100311, 100321,
100311, 100321, 100311, 100321, 100331, 100331,9020221,9020221,
10221, 10221,9030221,9030221};
_outgoingP = { -211, 111, 211, 311, 321, 321, 311,
111, 111, -211, 211, -321, -311, -321, -311,
-211, 111, -211, 111};
_outgoingV = { 213, 213, 113, 113, 113, -213, 213,
313, 323, 323, 313, 323, 313, 323, 313,
20213, 20113, 20213, 20113};
_coupling = { 3.57, 3.57, 3.57, 1., 1., 1.41, 1.41,
1.55, 1.55, 2.19, 2.19, 2.92, 2.92, 0.956, 0.956,
2.68, 2.68, 1.147, 1.147};
_maxweight = { 4.5, 4.5, 4.5, 4., 4., 4., 4.,
2., 2., 2., 2., 3.5, 3.5, 4., 4.,
4.5, 4.5, 3.2, 3.2};
// initial size of the arrays
_initsize=_incoming.size();
// intermediates
generateIntermediates(false);
}
void PScalarPScalarVectorDecayer::doinit() {
DecayIntegrator::doinit();
// check the parameters arew consistent
unsigned int isize=_coupling.size();
if(isize!=_incoming.size() || isize!=_outgoingP.size()||
isize!=_outgoingV.size() || isize!=_maxweight.size())
throw InitException() << "Inconsistent parameters in PScalarPScalarVectorDecayer"
<< Exception::abortnow;
// set up the integration channels
PhaseSpaceModePtr mode;
for(unsigned int ix=0;ix<_incoming.size();++ix) {
tPDPtr in = getParticleData( _incoming[ix]);
tPDVector out = {getParticleData(_outgoingP[ix]),
getParticleData(_outgoingV[ix])};
if(in&&out[0]&&out[1])
mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix]));
else
mode=PhaseSpaceModePtr();
addMode(mode);
}
}
int PScalarPScalarVectorDecayer::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 ==_outgoingP[ix]&&id2 ==_outgoingV[ix])||
(id2 ==_outgoingP[ix]&&id1 ==_outgoingV[ix])) imode=ix;
}
if(idbar==_incoming[ix]) {
if((id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix])||
(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix])) {
imode=ix;
cc=true;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
return imode;
}
void PScalarPScalarVectorDecayer::persistentOutput(PersistentOStream & os) const {
os << _coupling << _incoming << _outgoingP << _outgoingV << _maxweight;
}
void PScalarPScalarVectorDecayer::persistentInput(PersistentIStream & is, int) {
is >> _coupling >> _incoming >> _outgoingP >> _outgoingV >> _maxweight;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PScalarPScalarVectorDecayer,DecayIntegrator>
describeHerwigPScalarPScalarVectorDecayer("Herwig::PScalarPScalarVectorDecayer", "HwSMDecay.so");
void PScalarPScalarVectorDecayer::Init() {
static ClassDocumentation<PScalarPScalarVectorDecayer> documentation
("The PScalarPScalarVectorDecayer class is designed for"
" the decay of a pseduoscalar meson to two spin-1 particles.");
static ParVector<PScalarPScalarVectorDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&PScalarPScalarVectorDecayer::_incoming,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarPScalarVectorDecayer,int> interfaceOutgoingScalar
("OutgoingPScalar",
"The PDG code for the outgoing pseudoscalar meson",
&PScalarPScalarVectorDecayer::_outgoingP,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarPScalarVectorDecayer,int> interfaceOutgoingVector
("OutgoingVector",
"The PDG code for the outgoing vector meson",
&PScalarPScalarVectorDecayer::_outgoingV,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarPScalarVectorDecayer,double> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&PScalarPScalarVectorDecayer::_coupling,
0, 0, 0, 0., 100., false, false, true);
static ParVector<PScalarPScalarVectorDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&PScalarPScalarVectorDecayer::_maxweight,
0, 0, 0, 0., 100., false, false, true);
}
void PScalarPScalarVectorDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
// set up the spin information for the decay products
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
ScalarWaveFunction::constructSpinInfo(decay[0],outgoing,true);
VectorWaveFunction::constructSpinInfo(_vectors,decay[1],
outgoing,true,false);
}
double PScalarPScalarVectorDecayer::me2(const int,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1)));
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
}
_vectors.resize(3);
bool massless = outgoing[1]->id()==ParticleID::gamma;
for(unsigned int ix=0;ix<3;++ix) {
if(massless && ix==1) continue;
_vectors[ix] = HelicityFunctions::polarizationVector(-momenta[1],ix,Helicity::outgoing);
}
// calculate the matrix element
Lorentz5Momentum psum(part.momentum()+momenta[0]);
for(unsigned int ix=0;ix<3;++ix) {
- (*ME())(0,0,ix)=_coupling[imode()]/part.mass()*(_vectors[ix]*psum);
+ (*ME())(0,0,ix) = Complex(_coupling[imode()]/part.mass()*(_vectors[ix]*psum));
}
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 = 4.*sqr(_coupling[imode()]*pcm/momenta[1].mass());
// cerr << "testing matrix element for " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << me << " " << (me-test)/(me+test) << "\n";
// output the answer
return me;
}
// specify the 1-2 matrix element to be used in the running width calculation
bool PScalarPScalarVectorDecayer::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(true);
do {
if(id ==_incoming[ix]) {
if(id1==_outgoingP[ix]&&id2==_outgoingV[ix]) {
imode=ix;
order=true;
}
if(id2==_outgoingP[ix]&&id1==_outgoingV[ix]) {
imode=ix;
order=false;
}
}
if(idbar==_incoming[ix]&&imode<0) {
if(id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix]) {
imode=ix;
order=true;
}
if(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix]) {
imode=ix;
order=false;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
coupling=_coupling[imode];
mecode=10;
return order;
}
// output the setup information for the particle database
void PScalarPScalarVectorDecayer::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() << ":OutgoingPScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "newdef " << name() << ":OutgoingVector " << ix << " "
<< _outgoingV[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix] << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "insert " << name() << ":OutgoingPScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "insert " << name() << ":OutgoingVector " << ix << " "
<< _outgoingV[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix] << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc b/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc
--- a/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc
+++ b/Decay/ScalarMeson/PScalarVectorVectorDecayer.cc
@@ -1,275 +1,275 @@
// -*- C++ -*-
//
// PScalarVectorVectorDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 PScalarVectorVectorDecayer class.
//
#include "PScalarVectorVectorDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void PScalarVectorVectorDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<_incoming.size();++ix)
_maxweight[ix] = mode(ix)->maxWeight();
}
}
PScalarVectorVectorDecayer::PScalarVectorVectorDecayer()
: _incoming(10), _outgoing1(10), _outgoing2(10), _coupling(10), _maxweight(10) {
// decay eta -> omega gamma
_incoming[0] = 331; _outgoing1[0] = 223; _outgoing2[0] = 22;
_coupling[0] = 0.1412/GeV; _maxweight[0] = 1.2;
// decay pi -> gamma gamma
_incoming[1] = 111; _outgoing1[1] = 22; _outgoing2[1] = 22;
_coupling[1] = 0.0178/GeV; _maxweight[1] = 1.1;
// decay eta -> gamma gamma
_incoming[2] = 221; _outgoing1[2] = 22; _outgoing2[2] = 22;
_coupling[2] = 0.0176/GeV; _maxweight[2] = 1.1;
// decay eta' -> gamma gamma
_incoming[3] = 331; _outgoing1[3] = 22; _outgoing2[3] = 22;
_coupling[3] = 0.0221/GeV; _maxweight[3] = 1.1;
// decay eta_c -> rho rho
_incoming[4] = 441; _outgoing1[4] = 213; _outgoing2[4] = -213;
_coupling[4] = 0.0525/GeV; _maxweight[4] = 2.7;
_incoming[5] = 441; _outgoing1[5] = 113; _outgoing2[5] = 113;
_coupling[5] = 0.0371/GeV; _maxweight[5] = 2.7;
// decay eta-c -> phi phi
_incoming[6] = 441; _outgoing1[6] = 333; _outgoing2[6] = 333;
_coupling[6] = 0.0267/GeV; _maxweight[6] = 9.;
// decay eta-c -> gamma gamma
_incoming[7] = 441; _outgoing1[7] = 22; _outgoing2[7] = 22;
_coupling[7] = 0.00521/GeV; _maxweight[7] = 1.2;
// decay eta_c -> K* K*
_incoming[8] = 441; _outgoing1[8] = 323; _outgoing2[8] = -323;
_coupling[8] = 0.0308/GeV; _maxweight[8] = 5.3;
_incoming[9] = 441; _outgoing1[9] = 313; _outgoing2[9] = -313;
_coupling[9] = 0.0308/GeV; _maxweight[9] = 5.3;
// initial size of the vectors
_initsize = _incoming.size();
// intermediates
generateIntermediates(false);
}
void PScalarVectorVectorDecayer::doinit() {
DecayIntegrator::doinit();
// check the parameters arew consistent
unsigned int isize(_coupling.size());
if(isize!=_incoming.size() || isize!=_outgoing1.size()||
isize!=_outgoing2.size() || isize!=_maxweight.size())
throw InitException() << "Inconsistent parameters in PScalarVectorVectorDecayer"
<< Exception::abortnow;
// set up the integration channels
vector<double> wgt;
PhaseSpaceModePtr mode;
for(unsigned int ix=0;ix<_incoming.size();++ix) {
tPDPtr in = getParticleData( _incoming[ix]);
tPDVector out = {getParticleData(_outgoing1[ix]),
getParticleData(_outgoing2[ix])};
if(in&&out[0]&&out[1])
mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix]));
else
mode=PhaseSpaceModePtr();
addMode(mode);
}
}
int PScalarVectorVectorDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
cc = false;
if(children.size()!=2) return -1;
int id(parent->id());
int id1(children[0]->id());
int id2(children[1]->id());
unsigned int ix(0);
int imode(-1);
do {
if(_incoming[ix]==id) {
if((id1==_outgoing1[ix]&&id2==_outgoing2[ix])||
(id2==_outgoing1[ix]&&id1==_outgoing2[ix])) imode=ix;
}
++ix;
}
while(imode<0&&ix<_incoming.size());
return imode;
}
void PScalarVectorVectorDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_coupling,1/GeV) << _incoming << _outgoing1 << _outgoing2 << _maxweight;
}
void PScalarVectorVectorDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_coupling,1/GeV) >> _incoming >> _outgoing1 >> _outgoing2 >> _maxweight;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PScalarVectorVectorDecayer,DecayIntegrator>
describeHerwigPScalarVectorVectorDecayer("Herwig::PScalarVectorVectorDecayer", "HwSMDecay.so");
void PScalarVectorVectorDecayer::Init() {
static ClassDocumentation<PScalarVectorVectorDecayer> documentation
("The PScalarVectorVectorDecayer class is designed for"
" the decay of a pseduoscalar meson to two spin-1 particles.");
static ParVector<PScalarVectorVectorDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&PScalarVectorVectorDecayer::_incoming,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarVectorVectorDecayer,int> interfaceOutcoming1
("FirstOutgoing",
"The PDG code for the first outgoing particle",
&PScalarVectorVectorDecayer::_outgoing1,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarVectorVectorDecayer,int> interfaceOutcoming2
("SecondOutgoing",
"The PDG code for the second outgoing particle",
&PScalarVectorVectorDecayer::_outgoing2,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<PScalarVectorVectorDecayer,InvEnergy> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&PScalarVectorVectorDecayer::_coupling,
1/GeV, 0, ZERO, ZERO, 10000/GeV, false, false, true);
static ParVector<PScalarVectorVectorDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&PScalarVectorVectorDecayer::_maxweight,
0, 0, 0, 0., 200., false, false, true);
}
void PScalarVectorVectorDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
// set up the spin information for the decay products
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
bool photon[2]={false,false};
for(unsigned int ix=0;ix<2;++ix)
photon[ix] = decay[ix]->id()==ParticleID::gamma;
for(unsigned int ix=0;ix<2;++ix)
VectorWaveFunction::constructSpinInfo(_vectors[ix],decay[ix],
outgoing,true,photon[ix]);
}
double PScalarVectorVectorDecayer::me2(const int,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin1)));
bool photon[2]={false,false};
for(unsigned int ix=0;ix<2;++ix)
photon[ix] = outgoing[ix]->id()==ParticleID::gamma;
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
}
for(unsigned int ix=0;ix<2;++ix) {
_vectors[ix].resize(3);
for(unsigned int ihel=0;ihel<3;++ihel) {
if(photon[ix] && ihel==1) continue;
_vectors[ix][ihel] = HelicityFunctions::polarizationVector(-momenta[ix],ihel,Helicity::outgoing);
}
}
// now compute the matrix element
InvEnergy2 fact(_coupling[imode()]/part.mass());
for(unsigned int ix=0;ix<3;++ix) {
if(photon[0] && ix==1) continue;
for(unsigned int iy=0;iy<3;++iy) {
if(photon[1] && iy==1) continue;
- (*ME())(0,ix,iy)=fact*epsilon(_vectors[0][ix],momenta[1],
- _vectors[1][iy])*momenta[0];
+ (*ME())(0,ix,iy)=Complex(fact*epsilon(_vectors[0][ix],momenta[1],
+ _vectors[1][iy])*momenta[0]);
}
}
double output = ME()->contract(_rho).real();
// test of the matrix element
// double test = 2.*sqr(fact*part.mass())*
// sqr(Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(),momenta[1].mass()));
// cerr << "testing the matrix element for " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << output << " " << (output-test)/(output+test) << "\n";
return output;
}
// specify the 1-2 matrix element to be used in the running width calculation
bool PScalarVectorVectorDecayer::twoBodyMEcode(const DecayMode & dm, int & itype,
double & coupling) const {
int imode(-1);
int id(dm.parent()->id());
ParticleMSet::const_iterator pit(dm.products().begin());
int id1((**pit).id());++pit;
int id2((**pit).id());
unsigned int ix(0);
do {
if(_incoming[ix]==id) {
if((id1==_outgoing1[ix]&&id2==_outgoing2[ix])||
(id2==_outgoing1[ix]&&id1==_outgoing2[ix])) imode=ix;
}
++ix;
}
while(imode<0&&ix<_incoming.size());
coupling=_coupling[imode]*dm.parent()->mass();
itype = 3;
return id1==_outgoing1[imode]&&id2==_outgoing2[imode];
}
// output the setup info for the particle database
void PScalarVectorVectorDecayer::dataBaseOutput(ofstream & output,
bool header) const {
if(header) output << "update decayers set parameters=\"";
// parameters for the DecayIntegrator base class
DecayIntegrator::dataBaseOutput(output,false);
for(unsigned int ix=0;ix<_incoming.size();++ix) {
if(ix<_initsize) {
output << "newdef " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "newdef " << name() << ":FirstOutgoing " << ix << " "
<< _outgoing1[ix] << "\n";
output << "newdef " << name() << ":SecondOutgoing " << ix << " "
<< _outgoing2[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "insert " << name() << ":FirstOutgoing " << ix << " "
<< _outgoing1[ix] << "\n";
output << "insert " << name() << ":SecondOutgoing " << ix << " "
<< _outgoing2[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc b/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc
--- a/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc
+++ b/Decay/ScalarMeson/ScalarMesonTensorScalarDecayer.cc
@@ -1,281 +1,281 @@
// -*- C++ -*-
//
// ScalarMesonTensorScalarDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 ScalarMesonTensorScalarDecayer class.
//
#include "ScalarMesonTensorScalarDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void ScalarMesonTensorScalarDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<_incoming.size();++ix)
if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight();
}
}
ScalarMesonTensorScalarDecayer::ScalarMesonTensorScalarDecayer()
: _incoming(3), _outgoingT(3), _outgoingS(3), _coupling(3), _maxweight(3) {
// D+ -> f_2 pi
_incoming[0] = 411; _outgoingT[0] = 225; _outgoingS[0] = 211;
_coupling[0] = 8.23E-7/GeV; _maxweight[0] = 5;
// chi_c0 -> K*_0 K*_2
_incoming[1] = 10441; _outgoingT[1] = 325; _outgoingS[1] = -10321;
_coupling[1] = 0.0217/GeV; _maxweight[1] = 5;
_incoming[2] = 10441; _outgoingT[2] = 315; _outgoingS[2] = -10311;
_coupling[2] = 0.0217/GeV; _maxweight[2] = 5;
// initial size of the arrays
_initsize=_incoming.size();
// intermediates
generateIntermediates(false);
}
void ScalarMesonTensorScalarDecayer::doinit() {
DecayIntegrator::doinit();
// check the parameters arew consistent
unsigned int isize=_coupling.size();
if(isize!=_incoming.size() || isize!=_outgoingT.size()||
isize!=_outgoingS.size() || isize!=_maxweight.size())
throw InitException() << "Inconsistent parameters in "
<< "ScalarMesonTensorScalarDecayer"
<< Exception::abortnow;
// set up the integration channels
for(unsigned int ix=0;ix<_incoming.size();++ix) {
tPDPtr in = getParticleData(_incoming[ix]);
tPDVector out = {getParticleData(_outgoingT[ix]),
getParticleData(_outgoingS[ix])};
PhaseSpaceModePtr mode;
if(in&&out[0]&&out[1])
mode = new_ptr(PhaseSpaceMode(in,out,_maxweight[ix]));
else
mode=PhaseSpaceModePtr();
addMode(mode);
}
}
int ScalarMesonTensorScalarDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
if(children.size()!=2) return -1;
int id0(parent->id());
int id0bar = parent->CC() ? parent->CC()->id() : id0;
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;
unsigned int ix(0);
int imode(-1);
do {
if(id0 ==_incoming[ix]) {
if((id1 ==_outgoingT[ix]&&id2 ==_outgoingS[ix])||
(id2 ==_outgoingT[ix]&&id1 ==_outgoingS[ix])) {
imode=ix;
cc=false;
}
}
if(id0bar==_incoming[ix]&&imode<0) {
if((id1bar==_outgoingT[ix]&&id2bar==_outgoingS[ix])||
(id2bar==_outgoingT[ix]&&id1bar==_outgoingS[ix])) {
imode=ix;
cc=true;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
return imode;
}
void ScalarMesonTensorScalarDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(_coupling,1/GeV) << _incoming << _outgoingT << _outgoingS << _maxweight;
}
void ScalarMesonTensorScalarDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(_coupling,1/GeV) >> _incoming >> _outgoingT >> _outgoingS >> _maxweight;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<ScalarMesonTensorScalarDecayer,DecayIntegrator>
describeHerwigScalarMesonTensorScalarDecayer("Herwig::ScalarMesonTensorScalarDecayer", "HwSMDecay.so");
void ScalarMesonTensorScalarDecayer::Init() {
static ClassDocumentation<ScalarMesonTensorScalarDecayer> documentation
("The ScalarMesonTensorScalarDecayer class is designed for"
" the decay of a pseduoscalar meson to two spin-1 particles.");
static ParVector<ScalarMesonTensorScalarDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&ScalarMesonTensorScalarDecayer::_incoming,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<ScalarMesonTensorScalarDecayer,int> interfaceOutcomingT
("OutgoingTensor",
"The PDG code for the outgoing tensor",
&ScalarMesonTensorScalarDecayer::_outgoingT,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<ScalarMesonTensorScalarDecayer,int> interfaceOutcomingS
("OutgoingScalar",
"The PDG code for the outgoing scalar",
&ScalarMesonTensorScalarDecayer::_outgoingS,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<ScalarMesonTensorScalarDecayer,InvEnergy> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&ScalarMesonTensorScalarDecayer::_coupling,
1/GeV, 0, ZERO, ZERO, 100./GeV, false, false, true);
static ParVector<ScalarMesonTensorScalarDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&ScalarMesonTensorScalarDecayer::_maxweight,
0, 0, 0, 0., 100., false, false, true);
}
void ScalarMesonTensorScalarDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
// set up the spin information for the decay products
ScalarWaveFunction::constructSpinInfo(const_ptr_cast<tPPtr>(&part),
incoming,true);
TensorWaveFunction::constructSpinInfo(_tensors,decay[0],
outgoing,true,false);
ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
}
double ScalarMesonTensorScalarDecayer::me2(const int,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin0,PDT::Spin2,PDT::Spin0)));
if(meopt==Initialize) {
ScalarWaveFunction::
calculateWaveFunctions(_rho,const_ptr_cast<tPPtr>(&part),incoming);
}
TensorWaveFunction twave(momenta[0],outgoing[0],Helicity::outgoing);
_tensors.resize(5);
for(unsigned int ihel=0;ihel<5;++ihel) {
twave.reset(ihel);
_tensors[ihel] = twave.wave();
}
// calculate the matrix element
InvEnergy2 fact(_coupling[imode()]/part.mass());
LorentzPolarizationVectorE vtemp;
for(unsigned int ix=0;ix<5;++ix) {
vtemp = _tensors[ix]*part.momentum();
- (*ME())(0,ix,0) = fact * momenta[1].dot(vtemp);
+ (*ME())(0,ix,0) = Complex(fact * momenta[1].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 = 2.*pow<4,1>(pcm)*sqr(_coupling[imode()]*part.mass())/
// 3./pow<4,1>(momenta[0].mass());
// cerr << "testing matrix element for " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << me << " " << (me-test)/(me+test) << "\n";
// output the answer
return me;
}
// specify the 1-2 matrix element to be used in the running width calculation
bool ScalarMesonTensorScalarDecayer::twoBodyMEcode(const DecayMode & dm, int & itype,
double & coupling) const {
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);
int imode(-1);
do {
if(id ==_incoming[ix]) {
if(id1==_outgoingT[ix]&&id2==_outgoingS[ix]) {
imode=ix;
order=true;
}
if(id2==_outgoingT[ix]&&id1==_outgoingS[ix]) {
imode=ix;
order=false;
}
}
if(idbar==_incoming[ix]&&imode<0) {
if(id1bar==_outgoingT[ix]&&id2bar==_outgoingS[ix]) {
imode=ix;
order=true;
}
if(id2bar==_outgoingT[ix]&&id1bar==_outgoingS[ix]) {
imode=ix;
order=false;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
coupling=_coupling[imode]*dm.parent()->mass();
itype = 11;
return order;
}
// output the setup information for the particle database
void ScalarMesonTensorScalarDecayer::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() << ":OutgoingTensor " << ix << " "
<< _outgoingT[ix] << "\n";
output << "newdef " << name() << ":OutgoingScalar " << ix << " "
<< _outgoingS[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV << "\n";
output << "newdef " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
else {
output << "insert " << name() << ":Incoming " << ix << " "
<< _incoming[ix] << "\n";
output << "insert " << name() << ":OutgoingTensor " << ix << " "
<< _outgoingT[ix] << "\n";
output << "insert " << name() << ":OutgoingScalar " << ix << " "
<< _outgoingS[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/Tau/TauDecayer.cc b/Decay/Tau/TauDecayer.cc
--- a/Decay/Tau/TauDecayer.cc
+++ b/Decay/Tau/TauDecayer.cc
@@ -1,360 +1,360 @@
// -*- C++ -*-
//
// TauDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 TauDecayer class.
//
// Author: Peter Richardson
//
#include "TauDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Helicity/VectorSpinInfo.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "Herwig/Decay/DecayVertex.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include "ThePEG/Helicity/FermionSpinInfo.h"
#include "ThePEG/StandardModel/StandardModelBase.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void TauDecayer::doinit() {
DecayIntegrator::doinit();
// make sure the current got initialised
current_->init();
// set up the phase-space channels
tPDPtr tau = getParticleData(ParticleID::tauminus);
tPDPtr nu = getParticleData(ParticleID::nu_tau);
Energy mtau(tau->mass());
vector<double> channelwgts;
modeMap_.clear();
for(unsigned int ix=0;ix<current_->numberOfModes();++ix) {
// get the external particles for this mode
tPDVector out = {nu};
int iq(0),ia(0);
tPDVector ptemp = current_->particles(-3,ix,iq,ia);
out.insert(std::end(out), std::begin(ptemp), std::end(ptemp));
// the maximum weight
double maxweight = wgtMax_.size()>numberModes() ? wgtMax_[numberModes()] : 1.;
// create the mode
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(tau,out,maxweight));
// create the first piece of the channel
PhaseSpaceChannel channel((PhaseSpaceChannel(mode),0,1));
if(!current_->createMode(-3,tcPDPtr(),FlavourInfo(),ix,mode,1,0,channel,mtau)) continue;
// the channel weights
// the weights for the channel
if(wgtLoc_.size()>numberModes()&&
wgtLoc_[numberModes()]+mode->channels().size()<=weights_.size()) {
channelwgts=vector<double>(weights_.begin()+wgtLoc_[numberModes()],
weights_.begin()+wgtLoc_[numberModes()]+mode->channels().size());
}
else {
channelwgts.resize(mode->channels().size(),1./(mode->channels().size()));
}
modeMap_.push_back(ix);
// special for the two body modes
if(out.size()==2) {
channelwgts.clear();
mode=new_ptr(PhaseSpaceMode(tau,out,maxweight));
}
mode->setWeights(channelwgts);
// need to do the weights
addMode(mode);
}
current_->reset();
current_->touch();
current_->update();
}
void TauDecayer::doinitrun() {
current_->initrun();
DecayIntegrator::doinitrun();
if(initialize()) {
weights_.clear();
wgtLoc_.clear();
wgtMax_.clear();
for(unsigned int ix=0;ix<numberModes();++ix) {
wgtMax_.push_back(mode(ix)->maxWeight());
wgtLoc_.push_back(weights_.size());
for(unsigned int iy=0;iy<mode(ix)->channels().size();++iy) {
weights_.push_back(mode(ix)->channels()[iy].weight());
}
}
}
}
bool TauDecayer::accept(tcPDPtr parent, const tPDVector & children) const {
bool allowed(false);
// find the neutrino
int idnu(0),idtemp,idin(parent->id());
vector<int> idother;
tPDVector::const_iterator pit = children.begin();
tPDVector::const_iterator pend = children.end();
for( ; pit!=pend;++pit) {
idtemp=(**pit).id();
if(abs(idtemp)==16) idnu=idtemp;
else idother.push_back(idtemp);
}
if((idnu==ParticleID::nu_tau && idin==ParticleID::tauminus)||
(idnu==ParticleID::nu_taubar && idin==ParticleID::tauplus )) {
allowed=current_->accept(idother);
}
return allowed;
}
int TauDecayer::modeNumber(bool & cc,tcPDPtr parent, const tPDVector & children) const {
int imode(-1);
tPDVector::const_iterator pit = children.begin();
tPDVector::const_iterator pend = children.end();
int idtemp;vector<int> idother;
for( ; pit!=pend;++pit) {
idtemp=(**pit).id();
if(abs(idtemp)!=16) idother.push_back(idtemp);
}
unsigned int itemp=current_->decayMode(idother);
for(unsigned int ix=0;ix<modeMap_.size();++ix) {
if(modeMap_[ix]==itemp) imode=ix;
}
// perform the decay
cc=parent->id()==ParticleID::tauplus;
return imode;
}
void TauDecayer::persistentOutput(PersistentOStream & os) const {
os << modeMap_ << current_ << wgtLoc_
<< wgtMax_ << weights_ << polOpt_ << tauMpol_ << tauPpol_;
}
void TauDecayer::persistentInput(PersistentIStream & is, int) {
is >> modeMap_ >> current_ >> wgtLoc_
>> wgtMax_ >> weights_ >> polOpt_ >> tauMpol_ >> tauPpol_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TauDecayer,DecayIntegrator>
describeHerwigTauDecayer("Herwig::TauDecayer", "HwTauDecay.so");
void TauDecayer::Init() {
static ClassDocumentation<TauDecayer> documentation
("The TauDecayer class is designed to use a weak current"
" to perform the decay of the tau.");
static Reference<TauDecayer,WeakCurrent> interfaceWeakCurrent
("WeakCurrent",
"The reference for the decay current to be used.",
&TauDecayer::current_, false, false, true, false, false);
static ParVector<TauDecayer,int> interfaceWeightLocation
("WeightLocation",
"The locations of the weights for a given channel in the vector",
&TauDecayer::wgtLoc_,
0, 0, 0, 0, 10000, false, false, true);
static ParVector<TauDecayer,double> interfaceWeightMax
("MaximumWeight",
"The maximum weight for a given channel.",
&TauDecayer::wgtMax_,
0, 0, 0, 0., 100., false, false, true);
static ParVector<TauDecayer,double> interfaceWeights
("Weights",
"The weights for the integration.",
&TauDecayer::weights_,
0, 0, 0, 0., 1., false, false, true);
static Switch<TauDecayer,bool> interfacePolarizationOption
("PolarizationOption",
"Option of forcing the polarization of the tau leptons, N.B. you"
" should only use this option for making distributions for"
" comparision if you really know what you are doing.",
&TauDecayer::polOpt_, false, false, false);
static SwitchOption interfacePolarizationOptionDefault
(interfacePolarizationOption,
"Default",
"Don't force the polarization use the full spin density matrices"
" to get the right answer",
false);
static SwitchOption interfacePolarizationOptionForce
(interfacePolarizationOption,
"Force",
"Force the polarizations",
true);
static Parameter<TauDecayer,double> interfaceTauMinusPolarization
("TauMinusPolarization",
"The polarization of the tau-, left=-1, right=+1 if this is forced.",
&TauDecayer::tauMpol_, 0.0, -1.0, 1.0,
false, false, Interface::limited);
static Parameter<TauDecayer,double> interfaceTauPlusPolarization
("TauPlusPolarization",
"The polarization of the tau+, left=-1, right=+1 if this is forced.",
&TauDecayer::tauPpol_, 0.0, -1.0, 1.0,
false, false, Interface::limited);
}
void TauDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
if(part.id()==ParticleID::tauminus) {
SpinorWaveFunction ::
constructSpinInfo(inSpin_,const_ptr_cast<tPPtr>(&part),incoming,true);
SpinorBarWaveFunction::
constructSpinInfo(inBar_,decay[0],outgoing,true);
}
else {
SpinorBarWaveFunction::
constructSpinInfo(inBar_ ,const_ptr_cast<tPPtr>(&part),incoming,true);
SpinorWaveFunction::
constructSpinInfo(inSpin_,decay[0],outgoing,true);
}
current_->constructSpinInfo(ParticleVector(decay.begin()+1,decay.end()));
}
// combine the currents to give the matrix element
double TauDecayer::me2(const int ichan, const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
// map the mode to those in the current
int mode(modeMap_[imode()]);
// extract info on the decaying particle
if(meopt==Initialize) {
// spin density matrix for the decaying particle
rho_ = RhoDMatrix(PDT::Spin1Half);
if(part.id()==ParticleID::tauminus)
SpinorWaveFunction ::calculateWaveFunctions(inSpin_,rho_,
const_ptr_cast<tPPtr>(&part),
incoming);
else
SpinorBarWaveFunction::calculateWaveFunctions(inBar_ ,rho_,
const_ptr_cast<tPPtr>(&part),
incoming);
// fix rho if no correlations
fixRho(rho_);
if(polOpt_) {
rho_(0,1) = rho_(1,0) = 0.;
if(part.id()==ParticleID::tauminus) {
rho_(0,0) = 0.5*(1.-tauMpol_);
rho_(1,1) = 0.5*(1.+tauMpol_);
}
else {
rho_(0,0) = 0.5*(1.+tauPpol_);
rho_(1,1) = 0.5*(1.-tauPpol_);
}
}
// work out the mapping for the hadron vector
constants_ = vector<unsigned int>(outgoing.size()+1);
iSpin_ = vector<PDT::Spin >(outgoing.size());
int itemp(1);
unsigned int ix(outgoing.size());
do {
--ix;
iSpin_[ix] = outgoing[ix]->iSpin();
itemp *= iSpin_[ix];
constants_[ix] = itemp;
}
while(ix>0);
constants_[outgoing.size()] = 1;
constants_[0 ] = constants_[1];
}
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,iSpin_)));
// calculate the spinors for the decay products
if(part.id()==ParticleID::tauminus) {
inBar_.resize(2);
for(unsigned int ihel=0;ihel<2;++ihel)
inBar_[ihel] = HelicityFunctions::dimensionedSpinorBar(-momenta[0],ihel,Helicity::outgoing);
}
else {
inSpin_.resize(2);
for(unsigned int ihel=0;ihel<2;++ihel)
inSpin_[ihel] = HelicityFunctions::dimensionedSpinor (-momenta[0],ihel,Helicity::outgoing);
}
// calculate the hadron current
Energy q;
vector<LorentzPolarizationVectorE>
hadron(current_->current(tcPDPtr(),FlavourInfo(),mode,ichan,q,tPDVector(outgoing.begin()+1,outgoing.end()),
vector<Lorentz5Momentum>(momenta.begin()+1,momenta.end()),meopt));
// prefactor
double pre = sqr(pow(part.mass()/q,int(outgoing.size()-3)));
// calculate the lepton current
LorentzPolarizationVectorE lepton[2][2];
for(unsigned ix=0;ix<2;++ix) {
for(unsigned iy=0;iy<2;++iy) {
if(part.id()==15)
lepton[ix][iy]=2.*inSpin_[ix].leftCurrent(inBar_[iy]);
else
lepton[iy][ix]=2.*inSpin_[ix].leftCurrent(inBar_[iy]);
}
}
// compute the matrix element
vector<unsigned int> ihel(outgoing.size()+1);
for(unsigned int hhel=0;hhel<hadron.size();++hhel) {
// map the index for the hadrons to a helicity state
for(unsigned int ix=outgoing.size();ix>1;--ix) {
ihel[ix]=(hhel%constants_[ix-1])/constants_[ix];
}
// loop over the helicities of the tau and neutrino and set up the matrix
// element
for(ihel[1]=0;ihel[1]<2;++ihel[1]){
for(ihel[0]=0;ihel[0]<2;++ihel[0]) {
- (*ME())(ihel)= lepton[ihel[0]][ihel[1]].dot(hadron[hhel])*
- SM().fermiConstant();
+ (*ME())(ihel)= Complex(lepton[ihel[0]][ihel[1]].dot(hadron[hhel])*
+ SM().fermiConstant());
}
}
}
// multiply by the CKM element
int iq,ia;
current_->decayModeInfo(mode,iq,ia);
double ckm(1.);
if(iq<=6) {
if(iq%2==0) ckm = SM().CKM(iq/2-1,(abs(ia)-1)/2);
else ckm = SM().CKM(abs(ia)/2-1,(iq-1)/2);
}
return 0.5*pre*ckm*(ME()->contract(rho_)).real();
}
// output the setup information for the particle database
void TauDecayer::dataBaseOutput(ofstream & output,bool header) const {
unsigned int ix;
if(header) output << "update decayers set parameters=\"";
DecayIntegrator::dataBaseOutput(output,false);
for(ix=0;ix<wgtLoc_.size();++ix) {
output << "insert " << name() << ":WeightLocation " << ix << " "
<< wgtLoc_[ix] << "\n";
}
for(ix=0;ix<wgtMax_.size();++ix) {
output << "insert " << name() << ":MaximumWeight " << ix << " "
<< wgtMax_[ix] << "\n";
}
for(ix=0;ix<weights_.size();++ix) {
output << "insert " << name() << ":Weights " << ix << " "
<< weights_[ix] << "\n";
}
current_->dataBaseOutput(output,false,true);
output << "newdef " << name() << ":WeakCurrent " << current_->name() << " \n";
output << "\n\" where BINARY ThePEGName=\"" << fullName() << "\";\n";
}
diff --git a/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc b/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc
--- a/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc
+++ b/Decay/TensorMeson/TensorMesonVectorPScalarDecayer.cc
@@ -1,356 +1,356 @@
// -*- C++ -*-
//
// TensorMesonVectorPScalarDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 TensorMesonVectorPScalarDecayer class.
//
#include "TensorMesonVectorPScalarDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/WaveFunction/TensorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/epsilon.h"
#include "Herwig/Decay/TwoBodyDecayMatrixElement.h"
#include "ThePEG/Helicity/HelicityFunctions.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void TensorMesonVectorPScalarDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()){
for(unsigned int ix=0;ix<_incoming.size();++ix)
if(mode(ix)) _maxweight[ix] = mode(ix)->maxWeight();
}
}
TensorMesonVectorPScalarDecayer::TensorMesonVectorPScalarDecayer()
: _incoming(31), _outgoingV(31), _outgoingP(31),
_coupling(31), _maxweight(31) {
// a_2 -> rho pi
_incoming[0] = 115; _outgoingV[0] = 213; _outgoingP[0] = -211;
_coupling[0] = 21.1/GeV2; _maxweight[0] = 10.;
_incoming[1] = 215; _outgoingV[1] = 113; _outgoingP[1] = 211;
_coupling[1] = 21.1/GeV2; _maxweight[1] = 9.;
_incoming[2] = 215; _outgoingV[2] = 213; _outgoingP[2] = 111;
_coupling[2] = 21.1/GeV2; _maxweight[2] = 9.;
// a_2+/- -> gamma pi+/-
_incoming[3] = 215; _outgoingV[3] = 22; _outgoingP[3] = 211;
_coupling[3] = 0.551/GeV2; _maxweight[3] = 2.;
// k_2 -> K_2 omega
_incoming[4] = 315; _outgoingV[4] = 223; _outgoingP[4] = 311;
_coupling[4] = 11.66/GeV2; _maxweight[4] = 17.;
_incoming[5] = 325; _outgoingV[5] = 223; _outgoingP[5] = 321;
_coupling[5] = 11.66/GeV2; _maxweight[5] = 20.5;
// k_2+/- -> K+/- gamma
_incoming[6] = 325; _outgoingV[6] = 22; _outgoingP[6] = 321;
_coupling[6] = 0.553/GeV2; _maxweight[6] = 2.2;
// B_c2 -> B_c gamma
_incoming[7] = 545; _outgoingV[7] = 22; _outgoingP[7] = 541;
_coupling[7] = 0.651/GeV2; _maxweight[7] = 2.;
// K_2 -> K rho
_incoming[8] = 325; _outgoingV[8] = 113; _outgoingP[8] = 321;
_coupling[8] = 10.14/GeV2; _maxweight[8] = 9.;
_incoming[9] = 325; _outgoingV[9] = 213; _outgoingP[9] = 311;
_coupling[9] = 14.33/GeV2; _maxweight[9] = 9.;
_incoming[10] = 315; _outgoingV[10] = 113; _outgoingP[10] = 311;
_coupling[10] = 10.14/GeV2; _maxweight[10] = 9.;
_incoming[11] = 315; _outgoingV[11] = -213; _outgoingP[11] = 321;
_coupling[11] = 14.33/GeV2; _maxweight[11] = 9.;
// K_2 -> K* pi
_incoming[12] = 325; _outgoingV[12] = 323; _outgoingP[12] = 111;
_coupling[12] = 9.733/GeV2; _maxweight[12] = 13;
_incoming[13] = 325; _outgoingV[13] = 313; _outgoingP[13] = 211;
_coupling[13] = 13.77/GeV2; _maxweight[13] = 11;
_incoming[14] = 315; _outgoingV[14] = 313; _outgoingP[14] = 111;
_coupling[14] = 9.733/GeV2; _maxweight[14] = 8.;
_incoming[15] = 315; _outgoingV[15] = 323; _outgoingP[15] = -211;
_coupling[15] = 13.77/GeV2; _maxweight[15] = 8.;
// D_2 -> D* pi
_incoming[16] = 425; _outgoingV[16] = 423; _outgoingP[16] = 111;
_coupling[16] = 8.035/GeV2; _maxweight[16] = 2.2;
_incoming[17] = 425; _outgoingV[17] = 413; _outgoingP[17] = -211;
_coupling[17] = 11.670/GeV2; _maxweight[17] = 2.4;
_incoming[18] = 415; _outgoingV[18] = 413; _outgoingP[18] = 111;
_coupling[18] = 6.801/GeV2; _maxweight[18] = 2.4;
_incoming[19] = 415; _outgoingV[19] = 423; _outgoingP[19] = 211;
_coupling[19] = 9.527/GeV2; _maxweight[19] = 2.;
// D_s2 -> D* K
_incoming[20] = 435; _outgoingV[20] = 423; _outgoingP[20] = 321;
_coupling[20] = 13.10/GeV2; _maxweight[20] = 2.2;
_incoming[21] = 435; _outgoingV[21] = 413; _outgoingP[21] = 311;
_coupling[21] = 13.10/GeV2; _maxweight[21] = 2.5;
// B_2 -> B* pi
_incoming[22] = 525; _outgoingV[22] = 523; _outgoingP[22] = 111;
_coupling[22] = 4.99/GeV2; _maxweight[22] = 2.1;
_incoming[23] = 525; _outgoingV[23] = 513; _outgoingP[23] = 211;
_coupling[23] = 7.059/GeV2; _maxweight[23] = 2.1;
_incoming[24] = 515; _outgoingV[24] = 513; _outgoingP[24] = 111;
_coupling[24] = 4.99/GeV2; _maxweight[24] = 2.1;
_incoming[25] = 515; _outgoingV[25] = 523; _outgoingP[25] = -211;
_coupling[25] = 7.059/GeV2; _maxweight[25] = 2.1;
// D_s2
_incoming[26] = 435; _outgoingV[26] = 423; _outgoingP[26] = 321;
_coupling[26] = 13.09/GeV2; _maxweight[26] = 2.2;
_incoming[27] = 435; _outgoingV[27] = 413; _outgoingP[27] = 311;
_coupling[27] = 13.09/GeV2; _maxweight[27] = 2.5;
// B_s2
_incoming[28] = 535; _outgoingV[28] = 523; _outgoingP[28] = -321;
_coupling[28] = 7.29/GeV2; _maxweight[28] = 2.4;
_incoming[29] = 535; _outgoingV[29] = 513; _outgoingP[29] = -311;
_coupling[29] = 9.43/GeV2; _maxweight[29] = 2.1;
// upsilon_2(1d) to chi_b gamma
_incoming[30] = 20555; _outgoingV[30] = 22; _outgoingP[30] = 10551;
_coupling[30] = 1.11/GeV2; _maxweight[30] = 2.4;
// initial size of the arrays
_initsize=_incoming.size();
// intermediates
generateIntermediates(false);
}
void TensorMesonVectorPScalarDecayer::doinit() {
DecayIntegrator::doinit();
// check consistence 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 TensorMesonVectorPScalarDecayer"
<< Exception::abortnow;
// 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);
}
}
int TensorMesonVectorPScalarDecayer::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 ==_outgoingP[ix]&&id2 ==_outgoingV[ix])||
(id2 ==_outgoingP[ix]&&id1 ==_outgoingV[ix])) imode=ix;
}
if(idbar==_incoming[ix]) {
if((id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix])||
(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix])) {
imode=ix;
cc=true;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
return imode;
}
void TensorMesonVectorPScalarDecayer::persistentOutput(PersistentOStream & os) const {
os << _incoming << _outgoingV << _outgoingP << _maxweight << ounit(_coupling,1/GeV2);
}
void TensorMesonVectorPScalarDecayer::persistentInput(PersistentIStream & is, int) {
is >> _incoming >> _outgoingV >> _outgoingP >> _maxweight >> iunit(_coupling,1/GeV2);
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TensorMesonVectorPScalarDecayer,DecayIntegrator>
describeHerwigTensorMesonVectorPScalarDecayer("Herwig::TensorMesonVectorPScalarDecayer", "HwTMDecay.so");
void TensorMesonVectorPScalarDecayer::Init() {
static ClassDocumentation<TensorMesonVectorPScalarDecayer> documentation
("The TensorMesonVectorPScalarDecayer class implements the"
" decay of a tensor meson to a spin-1 particle and a pseduoscalar meson");
static ParVector<TensorMesonVectorPScalarDecayer,int> interfaceIncoming
("Incoming",
"The PDG code for the incoming particle",
&TensorMesonVectorPScalarDecayer::_incoming,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<TensorMesonVectorPScalarDecayer,int> interfaceOutcomingV
("OutgoingVector",
"The PDG code for the outgoing spin-1particle",
&TensorMesonVectorPScalarDecayer::_outgoingV,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<TensorMesonVectorPScalarDecayer,int> interfaceOutcomingP
("OutgoingScalar",
"The PDG code for the outgoing pseudoscalar meson",
&TensorMesonVectorPScalarDecayer::_outgoingP,
0, 0, 0, -10000000, 10000000, false, false, true);
static ParVector<TensorMesonVectorPScalarDecayer,InvEnergy2> interfaceCoupling
("Coupling",
"The coupling for the decay mode",
&TensorMesonVectorPScalarDecayer::_coupling,
1/GeV2, 0, ZERO, ZERO, 100./GeV2, false, false, true);
static ParVector<TensorMesonVectorPScalarDecayer,double> interfaceMaxWeight
("MaxWeight",
"The maximum weight for the decay mode",
&TensorMesonVectorPScalarDecayer::_maxweight,
0, 0, 0, 0., 1000., false, false, true);
}
void TensorMesonVectorPScalarDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
TensorWaveFunction::constructSpinInfo(_tensors,const_ptr_cast<tPPtr>(&part),
incoming,true,false);
// set up the spin information for the decay products
VectorWaveFunction::constructSpinInfo(_vectors,decay[0],outgoing,true,
decay[0]->id()==ParticleID::gamma);
ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
}
// matrix elememt for the process
double TensorMesonVectorPScalarDecayer::me2(const int,const Particle & part,
const tPDVector &,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(TwoBodyDecayMatrixElement(PDT::Spin2,PDT::Spin1,PDT::Spin0)));
// check for photons
bool photon(_outgoingV[imode()]==ParticleID::gamma);
// stuff for incoming particle
if(meopt==Initialize) {
_rho = RhoDMatrix(PDT::Spin2);
TensorWaveFunction::
calculateWaveFunctions(_tensors,_rho,const_ptr_cast<tPPtr>(&part),
incoming,false);
}
_vectors.resize(3);
for(unsigned int ix=0;ix<3;++ix) {
if(photon && ix==1) continue;
_vectors[ix] = HelicityFunctions::polarizationVector(-momenta[0],ix,Helicity::outgoing);
}
InvEnergy3 fact(_coupling[imode()]/part.mass());
// calculate the matrix element
for(unsigned int inhel=0;inhel<5;++inhel) {
for(unsigned int vhel=0;vhel<3;++vhel){
if(vhel==1&&photon) (*ME())(inhel,vhel,0)=0.;
else {
LorentzVector<complex<InvEnergy> > vtemp=
fact*epsilon(momenta[0],_vectors[vhel],momenta[1]);
- (*ME())(inhel,vhel,0)= (momenta[1]*_tensors[inhel]).dot(vtemp);
+ (*ME())(inhel,vhel,0) = Complex((momenta[1]*_tensors[inhel]).dot(vtemp));
}
}
}
double me = ME()->contract(_rho).real();
// test of the answer
// Energy pcm = Kinematics::pstarTwoBodyDecay(part.mass(),momenta[0].mass(),
// momenta[1].mass());
// double test = Energy4(pow<4,1>(2*pcm))*sqr( _coupling[imode()])/80.;
// cout << "testing matrix element for " << part.PDGName() << " -> "
// << outgoing[0]->PDGName() << " " << outgoing[1]->PDGName() << " "
// << me << " " << test << " " << (me-test)/(me+test) << endl;
// return the answer
return me;
}
bool TensorMesonVectorPScalarDecayer::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==_outgoingP[ix]&&id2==_outgoingV[ix]) {
imode=ix;
order=true;
}
if(id2==_outgoingP[ix]&&id1==_outgoingV[ix]) {
imode=ix;
order=false;
}
}
if(idbar==_incoming[ix]&&imode<0) {
if(id1bar==_outgoingP[ix]&&id2bar==_outgoingV[ix]) {
imode=ix;
order=true;
}
if(id2bar==_outgoingP[ix]&&id1bar==_outgoingV[ix]) {
imode=ix;
order=false;
}
}
++ix;
}
while(ix<_incoming.size()&&imode<0);
coupling=_coupling[imode]*sqr(dm.parent()->mass());
mecode=8;
return order;
}
void TensorMesonVectorPScalarDecayer::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() << ":OutgoingScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "newdef " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV2 << "\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() << ":OutgoingScalar " << ix << " "
<< _outgoingP[ix] << "\n";
output << "insert " << name() << ":Coupling " << ix << " "
<< _coupling[ix]*GeV2 << "\n";
output << "insert " << name() << ":MaxWeight " << ix << " "
<< _maxweight[ix] << "\n";
}
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h
--- a/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h
+++ b/Decay/VectorMeson/VectorMeson2BaryonsDecayer.h
@@ -1,230 +1,230 @@
// -*- C++ -*-
#ifndef Herwig_VectorMeson2BaryonsDecayer_H
#define Herwig_VectorMeson2BaryonsDecayer_H
//
// This is the declaration of the VectorMeson2BaryonsDecayer class.
//
#include "Herwig/Decay/DecayIntegrator.h"
#include "Herwig/Decay/PhaseSpaceMode.h"
#include "ThePEG/Helicity/LorentzPolarizationVector.h"
#include "ThePEG/Helicity/LorentzSpinorBar.h"
#include "ThePEG/Helicity/LorentzRSSpinorBar.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
*
* The <code>VectorMeson2BaryonsDecayer</code> class is designed for the decay
* of a vector meson to a baryon-antibaryon pair.
*
* In this case the matrix element is taken to have the form
* \f[\mathcal{M} = e_g \epsilon_\mu \bar{u}(p_f) \left[\gamma^\mu \right]u(p_{\bar{f}}).\f]
*
* The incoming vector mesons together with their decay products and the coupling
* \f$e_g\f$, \f$G_E\f$ and the phase \f$\phi_E\f$ can be specified using the interfaces for the class.
* The maximum weights
* for the decays can be calculated using the Initialize interface of the
* DecayIntegrator class or specified using the interface.
*
* The incoming and outgoing particles, couplings and maximum weights for
* many of the common \f$V\to f\bar{f}\f$ decays are specified in the default
* constructor.
*
* @see DecayIntegrator
* @see \ref VectorMeson2BaryonsDecayerInterfaces "The interfaces"
* defined for VectorMeson2BaryonsDecayer.
*
*/
class VectorMeson2BaryonsDecayer: public DecayIntegrator {
public:
/**
* The default constructor.
*/
VectorMeson2BaryonsDecayer();
/**
* 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 outgoing The particles produced in the decay
* @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
double me2(const int ichan,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const;
/**
* Construct the SpinInfos for the particles produced in the decay
*/
virtual void constructSpinInfo(const Particle & part,
ParticleVector outgoing) const;
/**
* 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);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
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:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- VectorMeson2BaryonsDecayer & operator=(const VectorMeson2BaryonsDecayer &);
+ VectorMeson2BaryonsDecayer & operator=(const VectorMeson2BaryonsDecayer &) = delete;
private:
/**
* \f$G_M\f$ coupling
*/
vector<double> gm_;
/**
* \f$G_E\f$ coupling
*/
vector<double> ge_;
/**
* Relative phase of \f$G_E\f$/\f$G_M\f$
*/
vector<double> phi_;
/**
* the PDG codes for the incoming particles
*/
vector<int> incoming_;
/**
* the PDG codes for the outgoing fermion
*/
vector<int> outgoingf_;
/**
* the PDG codes for the outgoing antifermion.
*/
vector<int> outgoinga_;
/**
* maximum weight for a decay
*/
vector<double> maxweight_;
/**
* Initial size of the vectors
*/
unsigned int initsize_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Polarization vectors for the decaying particle
*/
mutable vector<Helicity::LorentzPolarizationVector> vectors_;
/**
* Spinors for the decay products
*/
mutable vector<Helicity::LorentzSpinor <SqrtEnergy> > wave_;
/**
* barred spinors for the decay products
*/
mutable vector<Helicity::LorentzSpinorBar<SqrtEnergy> > wavebar_;
/**
* Spinors for the decay products
*/
mutable vector<Helicity::LorentzRSSpinor <SqrtEnergy> > wave2_;
/**
* barred spinors for the decay products
*/
mutable vector<Helicity::LorentzRSSpinorBar<SqrtEnergy> > wave2bar_;
};
}
#endif /* Herwig_VectorMeson2BaryonsDecayer_H */
diff --git a/Decay/VectorMeson/f1FourPiDecayer.cc b/Decay/VectorMeson/f1FourPiDecayer.cc
--- a/Decay/VectorMeson/f1FourPiDecayer.cc
+++ b/Decay/VectorMeson/f1FourPiDecayer.cc
@@ -1,313 +1,313 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the f1FourPiDecayer class.
//
#include "f1FourPiDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/epsilon.h"
using namespace Herwig;
f1FourPiDecayer::f1FourPiDecayer() :
gRhoPiPi_(6.), ga1RhoPi_(4.8*GeV), gf1a1Pi_(9.77/GeV), maxWeight_({1.,1.}) {
generateIntermediates(true);
}
// normally not implemented but do it here to get rid of the a_1 which
// can't be on-shell
ParticleVector f1FourPiDecayer::decay(const Particle & parent,
const tPDVector & children) const {
ParticleVector output = DecayIntegrator::decay(parent,children);
ParticleVector::iterator it =output.begin();
for(;it!=output.end();++it) {
long id = (**it).id();
if(id==20113 || id==20213 || id==-20213)
break;
}
if(it!=output.end()) {
PPtr a1 = *it;
output.erase(it);
for(PPtr child : a1->children()) {
output.push_back(child);
a1->abandonChild(child);
}
}
return output;
}
IBPtr f1FourPiDecayer::clone() const {
return new_ptr(*this);
}
IBPtr f1FourPiDecayer::fullclone() const {
return new_ptr(*this);
}
void f1FourPiDecayer::persistentOutput(PersistentOStream & os) const {
os << gRhoPiPi_ << ounit(ga1RhoPi_,GeV) << ounit(gf1a1Pi_,1./GeV)
<< ounit(ma1_,GeV) << ounit(ga1_,GeV)
<< ounit(mrho_,GeV) << ounit(grho_,GeV) << maxWeight_;
}
void f1FourPiDecayer::persistentInput(PersistentIStream & is, int) {
is >> gRhoPiPi_ >> iunit(ga1RhoPi_,GeV) >> iunit(gf1a1Pi_,1./GeV)
>> iunit(ma1_,GeV) >> iunit(ga1_,GeV)
>> iunit(mrho_,GeV) >> iunit(grho_,GeV) >> maxWeight_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<f1FourPiDecayer,DecayIntegrator>
describeHerwigf1FourPiDecayer("Herwig::f1FourPiDecayer", "HwVMDecay.so");
void f1FourPiDecayer::Init() {
static ClassDocumentation<f1FourPiDecayer> documentation
("The f1FourPiDecayer class implements a simple model for "
"f1 -> pipirho via an intermediate a_1");
static ParVector<f1FourPiDecayer,double> interfaceMaxWeight
("MaxWeight",
"Maximum weights for the decays",
&f1FourPiDecayer::maxWeight_, 2, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<f1FourPiDecayer,double> interfacegRhoPiPi
("gRhoPiPi",
"The coupling of the rho to two pions",
&f1FourPiDecayer::gRhoPiPi_, 6., 0.0, 10.0,
false, false, Interface::limited);
static Parameter<f1FourPiDecayer,Energy> interfacega1RhoPi
("ga1RhoPi",
"Coupling of the a_1 to rho pi",
&f1FourPiDecayer::ga1RhoPi_, GeV, 4.8*GeV, 0.*GeV, 20.*GeV,
false, false, Interface::limited);
static Parameter<f1FourPiDecayer,InvEnergy> interfacegf1a1Pi
("gf1a1Pi",
"Coupling of f_1 to a_1 pi",
&f1FourPiDecayer::gf1a1Pi_, 1./GeV, 1.0/GeV, 0.0/GeV, 10.0/GeV,
false, false, Interface::limited);
}
void f1FourPiDecayer::doinit() {
DecayIntegrator::doinit();
// pointers to the particles we need as external particles
tPDPtr f1 = getParticleData(ParticleID::f_1);
tPDPtr a1p = getParticleData(ParticleID::a_1plus);
tPDPtr a1m = getParticleData(ParticleID::a_1minus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
tPDPtr pip = getParticleData(ParticleID::piplus);
tPDPtr pim = getParticleData(ParticleID::piminus);
tPDPtr pi0 = getParticleData(ParticleID::pi0);
tPDPtr rhop = getParticleData(ParticleID::rhoplus);
tPDPtr rhom = getParticleData(ParticleID::rhominus);
tPDPtr rho0 = getParticleData(ParticleID::rho0);
// decay mode f_1 -> pi+ pi- pi+ pi-
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(f1,{pip,pim,pip,pim},maxWeight_[0]));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rho0,2,3,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,rho0,2,3,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,2,1,rho0,2,1,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,3,1,rho0,2,1,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rho0,2,3,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,1,1,rho0,2,3,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,4,1,rho0,2,1,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,3,1,rho0,2,1,2,2));
addMode(mode);
// decay mode f_1 -> pi+ pi0 pi- pi0
mode = new_ptr(PhaseSpaceMode(f1,{pip,pi0,pim,pi0},maxWeight_[0]));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rhop,2,3,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,rhop,2,3,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rhop,2,3,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,1,1,rhop,2,3,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,4,1,rhom,2,1,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,3,1,rhom,2,1,2,2));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,2,1,rhom,2,1,2,4));
mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,3,1,rhom,2,1,2,4));
addMode(mode);
// masses of intermediates
ma1_ = a1p->mass();
ga1_ = a1p->width();
mrho_ = rhop->mass();
grho_ = rhop->width();
}
void f1FourPiDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<maxWeight_.size();++ix)
maxWeight_[ix] = mode(ix)->maxWeight();
}
}
int f1FourPiDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
if(children.size()!=4 || parent->id()!=ParticleID::f_1) return -1;
// 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;
}
cc = false;
// f_1 -> 2pi+ 2pi-mode
if(npiplus==2 && npiminus==2) return 0;
// f_1 -> pi+ pi- pi0 pi0 mode
else if (npiplus ==1 && npiminus==1 && npi0==2) return 1;
// not found
return -1;
}
void f1FourPiDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(vector_,const_ptr_cast<tPPtr>(&part),
incoming,true,false);
// set up the spin information for the pions
for(unsigned int ix=0;ix<4;++ix)
ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
}
double f1FourPiDecayer::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,PDT::Spin0)));
useMe();
// polarization vectors for incoming
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vector_,rho_,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
// breit wigners
Lorentz5Momentum pa1[4] = {part.momentum()-momenta[0],part.momentum()-momenta[1],
part.momentum()-momenta[2],part.momentum()-momenta[3]};
complex<InvEnergy2> bwa1[4];
for(unsigned int ix=0;ix<4;++ix) {
pa1[ix].rescaleMass();
bwa1[ix] = Resonance::BreitWignera1(pa1[ix].mass2(),ma1_,ga1_)/sqr(ma1_);
}
// compute the matrix element (as a current and then contract)
LorentzVector<complex<InvEnergy> > current;
// decay mode f_1 -> pi+ pi- pi+pi-
double sym(0.5);
if(imode()==0) {
sym=0.25;
complex<InvEnergy2> bwrho[4] =
{Resonance::BreitWignerPWave((momenta[0]+momenta[1]).m2(),mrho_,grho_,momenta[0].mass(),momenta[1].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[0]+momenta[3]).m2(),mrho_,grho_,momenta[0].mass(),momenta[3].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[2]+momenta[1]).m2(),mrho_,grho_,momenta[2].mass(),momenta[1].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[2]+momenta[3]).m2(),mrho_,grho_,momenta[2].mass(),momenta[3].mass())/sqr(mrho_)};
if(ichan<=0) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rho0,2,3,2,4));
current += bwa1[0]*bwrho[3]*epsilon(part.momentum(),pa1[0],momenta[2]-momenta[3]);
}
if(ichan<0||ichan==1) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,rho0,2,3,2,4));
current -= bwa1[1]*bwrho[3]*epsilon(part.momentum(),pa1[1],momenta[2]-momenta[3]);
}
if(ichan<0||ichan==2) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,2,1,rho0,2,1,2,4));
current += bwa1[2]*bwrho[1]*epsilon(part.momentum(),pa1[2],momenta[0]-momenta[3]);
}
if(ichan<0||ichan==3) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,3,1,rho0,2,1,2,4));
current -= bwa1[1]*bwrho[1]*epsilon(part.momentum(),pa1[1],momenta[0]-momenta[3]);
}
if(ichan<0||ichan==4) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rho0,2,3,2,2));
current += bwa1[0]*bwrho[2]*epsilon(part.momentum(),pa1[0],momenta[2]-momenta[1]);
}
if(ichan<0||ichan==5) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,1,1,rho0,2,3,2,2));
current -= bwa1[3]*bwrho[2]*epsilon(part.momentum(),pa1[3],momenta[2]-momenta[1]);
}
if(ichan<0||ichan==6) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,3,1,4,1,rho0,2,1,2,2));
current += bwa1[2]*bwrho[0]*epsilon(part.momentum(),pa1[2],momenta[0]-momenta[1]);
}
if(ichan<0||ichan==7) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,4,1,3,1,rho0,2,1,2,2));
current -= bwa1[3]*bwrho[0]*epsilon(part.momentum(),pa1[3],momenta[0]-momenta[1]);
}
}
// decay mode f_1 -> pi+ pi0 pi- pi0
else {
complex<InvEnergy2> bwrho[4] =
{Resonance::BreitWignerPWave((momenta[0]+momenta[1]).m2(),mrho_,grho_,momenta[0].mass(),momenta[1].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[0]+momenta[3]).m2(),mrho_,grho_,momenta[0].mass(),momenta[3].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[2]+momenta[1]).m2(),mrho_,grho_,momenta[2].mass(),momenta[1].mass())/sqr(mrho_),
Resonance::BreitWignerPWave((momenta[2]+momenta[3]).m2(),mrho_,grho_,momenta[2].mass(),momenta[3].mass())/sqr(mrho_)};
double f1 = (momenta[2].mass2()-momenta[3].mass2())/sqr(mrho_);
double f2 = 1+f1;
f1 = 1.-f1;
if(ichan<=0) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,rhop,2,3,2,4));
current += bwa1[0]*bwrho[3]*epsilon(part.momentum(),pa1[0],f1*momenta[2]-f2*momenta[3]);
}
if(ichan<0||ichan==1) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,rhop,2,3,2,4));
current -= bwa1[1]*bwrho[3]*epsilon(part.momentum(),pa1[1],f1*momenta[2]-f2*momenta[3]);
}
if(ichan<0||ichan==2) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,4,1,rhop,2,3,2,2));
current += bwa1[0]*bwrho[2]*epsilon(part.momentum(),pa1[0],f1*momenta[2]-f2*momenta[1]);
}
if(ichan<0||ichan==3) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,1,1,rhop,2,3,2,2));
current -= bwa1[3]*bwrho[2]*epsilon(part.momentum(),pa1[3],f1*momenta[2]-f2*momenta[1]);
}
if(ichan<0||ichan==4) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,4,1,rhom,2,1,2,2));
current += bwa1[2]*bwrho[0]*epsilon(part.momentum(),pa1[2],f1*momenta[0]-f2*momenta[1]);
}
if(ichan<0||ichan==5) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,4,1,3,1,rhom,2,1,2,2));
current -= bwa1[3]*bwrho[0]*epsilon(part.momentum(),pa1[3],f1*momenta[0]-f2*momenta[1]);
}
if(ichan<0||ichan==6) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,3,1,2,1,rhom,2,1,2,4));
current += bwa1[2]*bwrho[1]*epsilon(part.momentum(),pa1[2],f1*momenta[0]-f2*momenta[3]);
}
if(ichan<0||ichan==7) {
// mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,3,1,rhom,2,1,2,4));
current -= bwa1[1]*bwrho[1]*epsilon(part.momentum(),pa1[1],f1*momenta[0]-f2*momenta[3]);
}
}
// contract the current
double pre = gRhoPiPi_*gf1a1Pi_*ga1RhoPi_;
for(unsigned int ihel=0;ihel<3;++ihel)
- (*ME())(ihel,0,0,0,0) = pre*current.dot(vector_[ihel])*part.mass();
+ (*ME())(ihel,0,0,0,0) = Complex(pre*current.dot(vector_[ihel])*part.mass());
// matrix element
return sym*ME()->contract(rho_).real();
}
void f1FourPiDecayer::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() << ":gRhoPiPi " << gRhoPiPi_ << "\n";
output << "newdef " << name() << ":ga1RhoPi " << ga1RhoPi_/GeV << "\n";
output << "newdef " << name() << ":gf1a1Pi " << gf1a1Pi_*GeV << "\n";
for(unsigned int ix=0;ix<maxWeight_.size();++ix) {
output << "newdef " << name() << ":maxWeight " << ix << " "
<< maxWeight_[ix] << "\n";
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/VectorMeson/f1RhoPiPiDecayer.cc b/Decay/VectorMeson/f1RhoPiPiDecayer.cc
--- a/Decay/VectorMeson/f1RhoPiPiDecayer.cc
+++ b/Decay/VectorMeson/f1RhoPiPiDecayer.cc
@@ -1,219 +1,219 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the f1RhoPiPiDecayer class.
//
#include "f1RhoPiPiDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/ParVector.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/epsilon.h"
using namespace Herwig;
f1RhoPiPiDecayer::f1RhoPiPiDecayer() :
ga1RhoPi_(4.8*GeV), gf1a1Pi_(9.77/GeV), maxWeight_({1.,1.}) {
generateIntermediates(true);
}
// normally not implemented but do it here to get rid of the a_1 which
// can't be on-shell
ParticleVector f1RhoPiPiDecayer::decay(const Particle & parent,
const tPDVector & children) const {
ParticleVector output = DecayIntegrator::decay(parent,children);
ParticleVector::iterator it =output.begin();
for(;it!=output.end();++it) {
long id = (**it).id();
if(id==20113 || id==20213 || id==-20213)
break;
}
if(it!=output.end()) {
PPtr a1 = *it;
output.erase(it);
for(PPtr child : a1->children()) {
output.push_back(child);
a1->abandonChild(child);
}
}
return output;
}
IBPtr f1RhoPiPiDecayer::clone() const {
return new_ptr(*this);
}
IBPtr f1RhoPiPiDecayer::fullclone() const {
return new_ptr(*this);
}
void f1RhoPiPiDecayer::persistentOutput(PersistentOStream & os) const {
os << ounit(ga1RhoPi_,GeV) << ounit(gf1a1Pi_,1./GeV)
<< ounit(ma1_,GeV) << ounit(ga1_,GeV) << maxWeight_;
}
void f1RhoPiPiDecayer::persistentInput(PersistentIStream & is, int) {
is >> iunit(ga1RhoPi_,GeV) >> iunit(gf1a1Pi_,1./GeV)
>> iunit(ma1_,GeV) >> iunit(ga1_,GeV) >> maxWeight_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<f1RhoPiPiDecayer,DecayIntegrator>
describeHerwigf1RhoPiPiDecayer("Herwig::f1RhoPiPiDecayer", "HwVMDecay.so");
void f1RhoPiPiDecayer::Init() {
static ClassDocumentation<f1RhoPiPiDecayer> documentation
("The f1RhoPiPiDecayer class implements a simple model for "
"f1 -> pipirho via an intermediate a_1");
static ParVector<f1RhoPiPiDecayer,double> interfaceMaxWeight
("MaxWeight",
"Maximum weights for the decays",
&f1RhoPiPiDecayer::maxWeight_, 2, 1.0, 0.0, 100.0,
false, false, Interface::limited);
static Parameter<f1RhoPiPiDecayer,Energy> interfacega1RhoPi
("ga1RhoPi",
"Coupling of the a_1 to rho pi",
&f1RhoPiPiDecayer::ga1RhoPi_, GeV, 4.8*GeV, 0.*GeV, 20.*GeV,
false, false, Interface::limited);
static Parameter<f1RhoPiPiDecayer,InvEnergy> interfacegf1a1Pi
("gf1a1Pi",
"Coupling of f_1 to a_1 pi",
&f1RhoPiPiDecayer::gf1a1Pi_, 1./GeV, 1.0/GeV, 0.0/GeV, 10.0/GeV,
false, false, Interface::limited);
}
void f1RhoPiPiDecayer::doinit() {
DecayIntegrator::doinit();
// pointers to the particles we need as external particles
tPDPtr f1 = getParticleData(ParticleID::f_1);
tPDPtr a1p = getParticleData(ParticleID::a_1plus);
tPDPtr a1m = getParticleData(ParticleID::a_1minus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
tPDPtr pip = getParticleData(ParticleID::piplus);
tPDPtr pim = getParticleData(ParticleID::piminus);
tPDPtr pi0 = getParticleData(ParticleID::pi0);
tPDPtr rhop = getParticleData(ParticleID::rhoplus);
tPDPtr rhom = getParticleData(ParticleID::rhominus);
tPDPtr rho0 = getParticleData(ParticleID::rho0);
// decay mode f_1 -> pi+ pi- rho0
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(f1,{pip,pim,rho0},maxWeight_[0]));
mode->addChannel((PhaseSpaceChannel(mode),0,a1m,0,1,1,2,1,3));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,2,1,1,1,3));
addMode(mode);
// decay mode f_1 -> pi- pi0 rho+
mode = new_ptr(PhaseSpaceMode(f1,{pim,pi0,rhop},maxWeight_[1]));
mode->addChannel((PhaseSpaceChannel(mode),0,a1p,0,1,1,2,1,3));
mode->addChannel((PhaseSpaceChannel(mode),0,a10,0,2,1,1,1,3));
addMode(mode);
ma1_ = a1p->mass();
ga1_ = a1p->width();
}
void f1RhoPiPiDecayer::doinitrun() {
DecayIntegrator::doinitrun();
if(initialize()) {
for(unsigned int ix=0;ix<maxWeight_.size();++ix)
maxWeight_[ix] = mode(ix)->maxWeight();
}
}
int f1RhoPiPiDecayer::modeNumber(bool & cc,tcPDPtr parent,
const tPDVector & children) const {
if(children.size()!=3 || parent->id()!=ParticleID::f_1) return -1;
// check the pions
int npi0(0),npiplus(0),npiminus(0),idrho(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;
else idrho=idtemp;
}
cc = false;
// f_1 -> pi+pi-rho0 mode
if(idrho==113 && npiplus==1 && npiminus==1) return 0;
// f_1 -> pi-pi0rho+ mode
else if (idrho== 213 && npiminus==1 && npi0==1) return 1;
else if (idrho==-213 && npiplus ==1 && npi0==1) {
cc=true;
return 1;
}
// not found
return -1;
}
void f1RhoPiPiDecayer::
constructSpinInfo(const Particle & part, ParticleVector decay) const {
VectorWaveFunction::constructSpinInfo(vectors_[0],const_ptr_cast<tPPtr>(&part),
incoming,true,false);
// set up the spin information for the decay products
// pions
for(unsigned int ix=0;ix<2;++ix)
ScalarWaveFunction::constructSpinInfo(decay[ix],outgoing,true);
// rho
VectorWaveFunction::constructSpinInfo(vectors_[1],decay[2],outgoing,true,false);
}
double f1RhoPiPiDecayer::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::Spin1)));
useMe();
// polarization vectors for incoming
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vectors_[0],rho_,
const_ptr_cast<tPPtr>(&part),
incoming,false);
}
// polarization vectors for rho
vectors_[1].resize(3);
for(unsigned int ix=0;ix<3;++ix)
vectors_[1][ix] = HelicityFunctions::polarizationVector(-momenta[2],ix,Helicity::outgoing);
// breit wigners
Lorentz5Momentum pa1[2] = {part.momentum()-momenta[0],
part.momentum()-momenta[1]};
for(unsigned int ix=0;ix<2;++ix) pa1[ix].rescaleMass();
complex<InvEnergy2> bwa1[2] = {Resonance::BreitWignera1(pa1[0].mass2(),ma1_,ga1_)/sqr(ma1_),
Resonance::BreitWignera1(pa1[1].mass2(),ma1_,ga1_)/sqr(ma1_)};
if(ichan>0) bwa1[ ichan == 0 ? 1 : 0 ] = ZERO;
// compute the matrix element
for(unsigned int ihel=0;ihel<3;++ihel) {
LorentzVector<complex<Energy2> > pol[2] = {epsilon(vectors_[0][ihel],part.momentum(),pa1[0]),
epsilon(vectors_[0][ihel],part.momentum(),pa1[1])};
for(unsigned int ohel=0;ohel<3;++ohel) {
- (*ME())(ihel,0,0,ohel) = gf1a1Pi_*ga1RhoPi_*(LorentzPolarizationVector(pol[0]*bwa1[0])-
- LorentzPolarizationVector(pol[1]*bwa1[1])).dot(vectors_[1][ohel]);
+ (*ME())(ihel,0,0,ohel) = Complex(gf1a1Pi_*ga1RhoPi_*(LorentzPolarizationVector(pol[0]*bwa1[0])-
+ LorentzPolarizationVector(pol[1]*bwa1[1])).dot(vectors_[1][ohel]));
}
}
// matrix element
return ME()->contract(rho_).real();
}
void f1RhoPiPiDecayer::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() << ":ga1RhoPi " << ga1RhoPi_/GeV << "\n";
output << "newdef " << name() << ":gf1a1Pi " << gf1a1Pi_*GeV << "\n";
for(unsigned int ix=0;ix<maxWeight_.size();++ix) {
output << "newdef " << name() << ":maxWeight " << ix << " "
<< maxWeight_[ix] << "\n";
}
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
diff --git a/Decay/WeakCurrents/ThreePionCLEOCurrent.cc b/Decay/WeakCurrents/ThreePionCLEOCurrent.cc
--- a/Decay/WeakCurrents/ThreePionCLEOCurrent.cc
+++ b/Decay/WeakCurrents/ThreePionCLEOCurrent.cc
@@ -1,1283 +1,1284 @@
// -*- C++ -*-
//
// ThreePionCLEOCurrent.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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 ThreePionCLEOCurrent class.
//
#include "ThreePionCLEOCurrent.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 "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/ResonanceHelpers.h"
#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
DescribeClass<ThreePionCLEOCurrent,WeakCurrent>
describeHerwigThreePionCLEOCurrent("Herwig::ThreePionCLEOCurrent",
"HwWeakCurrents.so");
HERWIG_INTERPOLATOR_CLASSDESC(ThreePionCLEOCurrent,Energy,Energy2)
ThreePionCLEOCurrent::ThreePionCLEOCurrent() {
addDecayMode(1,-1);
addDecayMode(2,-2);
addDecayMode(2,-1);
addDecayMode(1,-1);
addDecayMode(2,-2);
addDecayMode(2,-1);
setInitialModes(6);
// rho masses and widths
_rhomass = {0.7743*GeV,1.370*GeV};
_rhowidth = {0.1491*GeV,0.386*GeV};
// f_2 mass and width
_f2mass = 1.275*GeV;
_f2width = 0.185*GeV;
// f_0(1370) mass and width
_f0mass = 1.186*GeV;
_f0width = 0.350*GeV;
// sigma mass and width
_sigmamass = 0.860*GeV;
_sigmawidth = 0.880*GeV;
// a1 mass and width
_a1mass = 1.331*GeV;
_a1width = 0.814*GeV;
// parameters for the K K* contribution to the a_1 running width
_mKstar = 894*MeV;
_mK = 496*MeV;
_gammk = 3.32;
// pion decay constant
_fpi = 130.7*MeV/sqrt(2.);
// couplings and phases for the different channels
// p-wave rho and rho prime
using Constants::pi;
_rhomagP = {1.,0.12};
_rhophaseP = {0.,0.99*pi};
// d-wave rho and rho prime
_rhomagD = {0.37/GeV2,0.87/GeV2};
_rhophaseD = {-0.15*pi, 0.53*pi};
// f_2
_f2mag = 0.71/GeV2;
_f2phase = 0.56*pi;
_f2coup = ZERO;
// sigma
_sigmamag = 2.10;
_sigmaphase = 0.23*pi;
_sigmacoup = 0.;
// f_0
_f0mag = 0.77;
_f0phase = -0.54*pi;
_f0coup = 0.;
// initialize the a_1 width
_initializea1=false;
_a1opt=true;
double a1q2in[200]={0 ,15788.6,31577.3,47365.9,63154.6,78943.2,
94731.9,110521 ,126309 ,142098 ,157886 ,173675 ,
189464 ,205252 ,221041 ,236830 ,252618 ,268407 ,
284196 ,299984 ,315773 ,331562 ,347350 ,363139 ,
378927 ,394716 ,410505 ,426293 ,442082 ,457871 ,
473659 ,489448 ,505237 ,521025 ,536814 ,552603 ,
568391 ,584180 ,599969 ,615757 ,631546 ,647334 ,
663123 ,678912 ,694700 ,710489 ,726278 ,742066 ,
757855 ,773644 ,789432 ,805221 ,821010 ,836798 ,
852587 ,868375 ,884164 ,899953 ,915741 ,931530 ,
947319 ,963107 ,978896 ,994685 ,
1.01047e+06,1.02626e+06,1.04205e+06,1.05784e+06,
1.07363e+06,1.08942e+06,1.10521e+06,1.12099e+06,
1.13678e+06,1.15257e+06,1.16836e+06,1.18415e+06,
1.19994e+06,1.21573e+06,1.23151e+06,1.24730e+06,
1.26309e+06,1.27888e+06,1.29467e+06,1.31046e+06,
1.32625e+06,1.34203e+06,1.35782e+06,1.37361e+06,
1.38940e+06,1.40519e+06,1.42098e+06,1.43677e+06,
1.45256e+06,1.46834e+06,1.48413e+06,1.49992e+06,
1.51571e+06,1.53150e+06,1.54729e+06,1.56308e+06,
1.57886e+06,1.59465e+06,1.61044e+06,1.62623e+06,
1.64202e+06,1.65781e+06,1.67360e+06,1.68939e+06,
1.70517e+06,1.72096e+06,1.73675e+06,1.75254e+06,
1.76833e+06,1.78412e+06,1.79991e+06,1.81569e+06,
1.83148e+06,1.84727e+06,1.86306e+06,1.87885e+06,
1.89464e+06,1.91043e+06,1.92621e+06,1.94200e+06,
1.95779e+06,1.97358e+06,1.98937e+06,2.00516e+06,
2.02095e+06,2.03674e+06,2.05252e+06,2.06831e+06,
2.08410e+06,2.09989e+06,2.11568e+06,2.13147e+06,
2.14726e+06,2.16304e+06,2.17883e+06,2.19462e+06,
2.21041e+06,2.22620e+06,2.24199e+06,2.25778e+06,
2.27356e+06,2.28935e+06,2.30514e+06,2.32093e+06,
2.33672e+06,2.35251e+06,2.36830e+06,2.38409e+06,
2.39987e+06,2.41566e+06,2.43145e+06,2.44724e+06,
2.46303e+06,2.47882e+06,2.49461e+06,2.51039e+06,
2.52618e+06,2.54197e+06,2.55776e+06,2.57355e+06,
2.58934e+06,2.60513e+06,2.62092e+06,2.63670e+06,
2.65249e+06,2.66828e+06,2.68407e+06,2.69986e+06,
2.71565e+06,2.73144e+06,2.74722e+06,2.76301e+06,
2.77880e+06,2.79459e+06,2.81038e+06,2.82617e+06,
2.84196e+06,2.85774e+06,2.87353e+06,2.88932e+06,
2.90511e+06,2.92090e+06,2.93669e+06,2.95248e+06,
2.96827e+06,2.98405e+06,2.99984e+06,3.01563e+06,
3.03142e+06,3.04721e+06,3.06300e+06,3.07879e+06,
3.09457e+06,3.11036e+06,3.12615e+06,3.14194e+06};
double a1widthin[200]={0,0,0,0,0,0,0,0,
0,0,0,0.00021256,0.0107225,0.0554708,0.150142,0.303848,
0.522655,0.81121,1.1736,1.61381,2.13606,2.74499,3.44583,4.24454,
5.14795,6.16391,7.3014,8.57079,9.98398,11.5547,13.2987,15.2344,
17.3827,19.7683,22.4195,25.3695,28.6568,32.3264,36.4311,41.0322,
46.201,52.0203,58.5847,66.0011,74.3871,83.8666,94.5615,106.578,
119.989,134.807,150.968,168.315,186.615,205.576,224.893,244.28,
263.499,282.364,300.748,318.569,335.781,352.367,368.327,383.677,
398.438,412.638,426.306,439.472,452.167,464.421,476.263,487.719,
498.815,509.576,520.024,530.179,540.063,549.693,559.621,568.26,
577.229,586.005,594.604,603.035,611.314,619.447,627.446,635.321,
643.082,650.736,658.288,665.75,673.127,680.427,687.659,694.82,
701.926,708.977,715.983,722.944,729.862,736.752,743.619,750.452,
757.271,764.076,770.874,777.658,784.444,791.233,798.027,804.838,
811.649,818.485,825.342,832.224,839.139,846.082,853.059,860.079,
867.143,874.248,881.409,919.527,945.28,965.514,983.228,999.471,
1014.69,1029.15,1043.05,1056.49,1069.57,1082.36,1094.88,1107.2,
1120.89,1131.4,1143.33,1155.15,1166.92,1178.61,1190.27,1201.92,
1213.55,1225.18,1236.81,1250.06,1260.16,1271.86,1283.64,1295.46,
1307.36,1319.3,1331.34,1343.45,1355.64,1367.93,1380.31,1392.77,
1405.35,1418.03,1430.83,1443.75,1457.17,1469.94,1483.22,1496.64,
1510.18,1523.86,1537.67,1551.64,1565.72,1579.99,1594.38,1608.92,
1623.63,1642.08,1653.51,1668.69,1684.03,1699.53,1715.21,1731.04,
1747.05,1763.23,1779.59,1796.12,1812.83,1829.72,1846.79,1864.04,
1881.49,1899.11,1916.93,1934.93,1953.13,1971.52,1990.12,2008.89};
if(_a1runwidth.empty()) {
vector<double> tmp1(a1widthin,a1widthin+200);
std::transform(tmp1.begin(), tmp1.end(),
back_inserter(_a1runwidth),
[](double x){return x*GeV;});
vector<double> tmp2(a1q2in,a1q2in+200);
_a1runq2.clear();
std::transform(tmp2.begin(), tmp2.end(),
back_inserter(_a1runq2),
[](double x){return x*GeV2;});
}
// zero parameters which will be calculated later to avoid problems
_mpi0=ZERO;
_mpic=ZERO;
_fact=ZERO;
_maxmass=ZERO;
_maxcalc=ZERO;
}
void ThreePionCLEOCurrent::doinit() {
WeakCurrent::doinit();
// parameters for the breit-wigners
_mpic = getParticleData(ParticleID::piplus)->mass();
_mpi0 = getParticleData(ParticleID::pi0) ->mass();
// 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));
// overall coupling
_fact = 2.*sqrt(2.)/_fpi/3.;
// initialise the a_1 running width calculation
inita1Width(-1);
}
void ThreePionCLEOCurrent::persistentOutput(PersistentOStream & os) const {
os << ounit(_rhomass,GeV) << ounit(_rhowidth,GeV)
<< ounit(_f2mass,GeV) << ounit(_f2width,GeV)
<< ounit(_f0mass,GeV) << ounit(_f0width,GeV)
<< ounit(_sigmamass,GeV) << ounit(_sigmawidth,GeV)
<< ounit(_mpi0,GeV) << ounit(_mpic,GeV)
<< ounit(_fpi,GeV) << ounit(_fact,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
<< ounit(_a1mass,GeV) << ounit(_a1width,GeV) << ounit(_a1runwidth,GeV)
<< ounit(_a1runq2,GeV2) << _initializea1
<< ounit(_mKstar,GeV) << ounit(_mK,GeV) << _gammk << _a1opt
<< ounit(_maxmass,GeV) << ounit(_maxcalc,GeV) << _a1runinter;
}
void ThreePionCLEOCurrent::persistentInput(PersistentIStream & is, int) {
is >> iunit(_rhomass,GeV) >> iunit(_rhowidth,GeV)
>> iunit(_f2mass,GeV) >> iunit(_f2width,GeV)
>> iunit(_f0mass,GeV) >> iunit(_f0width,GeV)
>> iunit(_sigmamass,GeV) >> iunit(_sigmawidth,GeV)
>> iunit(_mpi0,GeV) >> iunit(_mpic,GeV)
>> iunit(_fpi,GeV) >> iunit(_fact,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
>> iunit(_a1mass,GeV) >> iunit(_a1width,GeV) >> iunit(_a1runwidth,GeV)
>> iunit(_a1runq2,GeV2) >> _initializea1
>> iunit(_mKstar,GeV) >> iunit(_mK,GeV) >> _gammk >> _a1opt
>> iunit(_maxmass,GeV) >> iunit(_maxcalc,GeV) >> _a1runinter;
}
void ThreePionCLEOCurrent::Init() {
static ClassDocumentation<ThreePionCLEOCurrent> documentation
("The ThreePionCLEOCurrent class performs the decay of the"
" tau to three pions using the currents from CLEO",
"The decay of tau to three pions is modelled using the currents from "
"\\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<ThreePionCLEOCurrent,Energy> interfacerhomass
("RhoMasses",
"The masses of the different rho resonnaces",
&ThreePionCLEOCurrent::_rhomass,
MeV, 0, ZERO, -10000*MeV, 10000*MeV, false, false, true);
static ParVector<ThreePionCLEOCurrent,Energy> interfacerhowidth
("RhoWidths",
"The widths of the different rho resonnaces",
&ThreePionCLEOCurrent::_rhowidth,
MeV, 0, ZERO, -10000*MeV, 10000*MeV, false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacef_2Mass
("f_2Mass",
"The mass of the f_2 meson",
&ThreePionCLEOCurrent::_f2mass, GeV, 1.275*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacef_2Width
("f_2Width",
"The width of the f_2 meson",
&ThreePionCLEOCurrent::_f2width, GeV, 0.185*GeV, ZERO, 1.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacef_0Mass
("f_0Mass",
"The mass of the f_0 meson",
&ThreePionCLEOCurrent::_f0mass, GeV, 1.186*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacef_0Width
("f_0Width",
"The width of the f_0 meson",
&ThreePionCLEOCurrent::_f0width, GeV, 0.350*GeV, ZERO, 1.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacesigmaMass
("sigmaMass",
"The mass of the sigma meson",
&ThreePionCLEOCurrent::_sigmamass, GeV, 0.860*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacesigmaWidth
("sigmaWidth",
"The width of the sigma meson",
&ThreePionCLEOCurrent::_sigmawidth, GeV, 0.880*GeV, ZERO, 2.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacea1Mass
("a1Mass",
"The mass of the a_1 meson",
&ThreePionCLEOCurrent::_a1mass, GeV, 1.331*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfacea1Width
("a1Width",
"The width of the a_1 meson",
&ThreePionCLEOCurrent::_a1width, GeV, 0.814*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfaceKaonMass
("KaonMass",
"The mass of the kaon",
&ThreePionCLEOCurrent::_mK, GeV, 0.496*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfaceKStarMass
("KStarMass",
"The mass of the k* meson",
&ThreePionCLEOCurrent::_mKstar, GeV, 0.894*GeV, ZERO, 10.0*GeV,
false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfaceKaonCoupling
("KaonCoupling",
"The relative coupling for the kaon in the a_1 running width",
&ThreePionCLEOCurrent::_gammk, 3.32, 0.0, 10.0,
false, false, true);
static Parameter<ThreePionCLEOCurrent,Energy> interfaceFpi
("Fpi",
"The pion decay constant",
&ThreePionCLEOCurrent::_fpi, MeV, 130.7*MeV/sqrt(2.), ZERO, 500.0*MeV,
false, false, true);
static ParVector<ThreePionCLEOCurrent,double> interfacerhomagP
("RhoPWaveMagnitude",
"The magnitude of the couplings for the p-wave rho currents",
&ThreePionCLEOCurrent::_rhomagP,
0, 0, 0, 0, 10000, false, false, true);
static ParVector<ThreePionCLEOCurrent,double> interfacerhophaseP
("RhoPWavePhase",
"The phase of the couplings for the p-wave rho currents",
&ThreePionCLEOCurrent::_rhophaseP,
0, 0, 0, -Constants::twopi, Constants::twopi, false, false, true);
static ParVector<ThreePionCLEOCurrent,InvEnergy2> interfacerhomagD
("RhoDWaveMagnitude",
"The magnitude of the couplings for the d-wave rho currents",
&ThreePionCLEOCurrent::_rhomagD,
1/MeV2, 0, ZERO, ZERO, 10000/MeV2, false, false, true);
static ParVector<ThreePionCLEOCurrent,double> interfacerhophaseD
("RhoDWavePhase",
"The phase of the couplings for the d-wave rho currents",
&ThreePionCLEOCurrent::_rhophaseD,
0, 0, 0, -Constants::twopi, Constants::twopi, false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfacef0Phase
("f0Phase",
"The phase of the f_0 scalar current",
&ThreePionCLEOCurrent::_f0phase, 0.54*Constants::pi, -Constants::twopi, Constants::twopi,
false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfacef2Phase
("f2Phase",
"The phase of the f_2 tensor current",
&ThreePionCLEOCurrent::_f2phase, 0.56*Constants::pi,-Constants::twopi, Constants::twopi,
false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfacesigmaPhase
("sigmaPhase",
"The phase of the sigma scalar current",
&ThreePionCLEOCurrent::_sigmaphase, 0.23*Constants::pi, -Constants::twopi, Constants::twopi,
false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfacef0Magnitude
("f0Magnitude",
"The magnitude of the f_0 scalar current",
&ThreePionCLEOCurrent::_f0mag, 0.77, 0.0, 10,
false, false, true);
static Parameter<ThreePionCLEOCurrent,InvEnergy2> interfacef2Magnitude
("f2Magnitude",
"The magnitude of the f_2 tensor current",
&ThreePionCLEOCurrent::_f2mag, 1./GeV2, 0.71/GeV2, ZERO, 10./GeV2,
false, false, true);
static Parameter<ThreePionCLEOCurrent,double> interfacesigmaMagnitude
("sigmaMagnitude",
"The magnitude of the sigma scalar current",
&ThreePionCLEOCurrent::_sigmamag, 2.1, 0.0, 10,
false, false, true);
static ParVector<ThreePionCLEOCurrent,Energy> interfacea1RunningWidth
("a1RunningWidth",
"The values of the a_1 width for interpolation to giving the running width.",
&ThreePionCLEOCurrent::_a1runwidth,
MeV, 0, ZERO, ZERO, 10000000*MeV, false, false, true);
static ParVector<ThreePionCLEOCurrent,Energy2> interfacea1RunningQ2
("a1RunningQ2",
"The values of the q^2 for interpolation to giving the running width.",
&ThreePionCLEOCurrent::_a1runq2,
MeV2, 0, ZERO, ZERO, 10000000*MeV2, false, false, true);
static Switch<ThreePionCLEOCurrent,bool> interfaceInitializea1
("Initializea1",
"Initialise the calculation of the a_1 running width",
&ThreePionCLEOCurrent::_initializea1, false, false, false);
static SwitchOption interfaceInitializea1Initialization
(interfaceInitializea1,
"Yes",
"Initialize the calculation",
true);
static SwitchOption interfaceInitializea1NoInitialization
(interfaceInitializea1,
"No",
"Use the default values",
false);
static Switch<ThreePionCLEOCurrent,bool> interfacea1WidthOption
("a1WidthOption",
"Option for the treatment of the a1 width",
&ThreePionCLEOCurrent::_a1opt, true, false, false);
static SwitchOption interfacea1WidthOptionLocal
(interfacea1WidthOption,
"Local",
"Use a calculation of the running width based on the parameters as"
" interpolation table.",
true);
static SwitchOption interfacea1WidthOptionParam
(interfacea1WidthOption,
"Kuhn",
"Use the parameterization of Kuhn and Santamaria for default parameters."
" This should only be used for testing vs TAUOLA",
false);
}
// initialisation of the a_1 width
// (iopt=-1 initialises, iopt=0 starts the interpolation)
void ThreePionCLEOCurrent::inita1Width(int iopt) {
if(iopt==-1) {
_maxcalc=_maxmass;
if(!_initializea1||_maxmass==ZERO) return;
// parameters for the table of values
Energy2 step=sqr(_maxmass)/200.;
// function to be integrated to give the matrix element
// 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);
Energy mrho=getParticleData(ParticleID::rhoplus)->mass();
Energy wrho=getParticleData(ParticleID::rhoplus)->width();
vector<Energy> inmass;inmass.push_back(mrho);inmass.push_back(mrho);
vector<Energy> inwidth;inwidth.push_back(wrho);inwidth.push_back(wrho);
vector<double> inpow(2,0.0);
ThreeBodyAllOnCalculator<ThreePionCLEOCurrent>
widthgenN(inweights,intype,inmass,inwidth,inpow,*this,0,_mpi0,_mpi0,_mpic);
ThreeBodyAllOnCalculator<ThreePionCLEOCurrent>
widthgenC(inweights,intype,inmass,inwidth,inpow,*this,1,_mpic,_mpic,_mpic);
// normalisation constant to give physical width if on shell
double a1const = _a1width/(widthgenN.partialWidth(sqr(_a1mass))+
widthgenC.partialWidth(sqr(_a1mass)));
// loop to give the values
_a1runq2.clear();_a1runwidth.clear();
for(Energy2 moff2=ZERO; moff2<=sqr(_maxmass); moff2+=step) {
Energy moff=sqrt(moff2);
_a1runq2.push_back(moff2);
Energy charged=a1const*widthgenC.partialWidth(moff2);
Energy neutral=a1const*widthgenN.partialWidth(moff2);
Energy kaon = moff<=_mK+_mKstar ? ZERO : 2.870*_gammk*_gammk/8./Constants::pi*
Kinematics::pstarTwoBodyDecay(moff,_mK,_mKstar)/moff2*GeV2;
Energy total = charged + neutral + kaon;
_a1runwidth.push_back(total);
}
}
// set up the interpolator
else if(iopt==0) {
_a1runinter = make_InterpolatorPtr(_a1runwidth,_a1runq2,3);
}
}
void ThreePionCLEOCurrent::CLEOFormFactor(int imode,int ichan,
Energy2 q2,Energy2 s1, Energy2 s2, Energy2 s3,
Complex & F1, Complex & F2,
Complex & F3) const {
useMe();
double fact=1.;
if(imode<=1) {
// identical particle factors
fact = 1./sqrt(6.);
// compute the breit wigners we need
Complex sigbws1 = Resonance::BreitWignerSWave(s1,_sigmamass,_sigmawidth,_mpi0,_mpi0);
Complex sigbws2 = Resonance::BreitWignerSWave(s2,_sigmamass,_sigmawidth,_mpi0,_mpi0);
Complex sigbws3 = Resonance::BreitWignerSWave(s3,_sigmamass,_sigmawidth,_mpi0,_mpi0);
Complex f0bws1 = Resonance::BreitWignerSWave(s1, _f0mass, _f0width,_mpi0,_mpi0);
Complex f0bws2 = Resonance::BreitWignerSWave(s2, _f0mass, _f0width,_mpi0,_mpi0);
Complex f0bws3 = Resonance::BreitWignerSWave(s3, _f0mass, _f0width,_mpi0,_mpi0);
Complex f2bws1 = Resonance::BreitWignerDWave(s1, _f2mass, _f2width,_mpi0,_mpi0);
Complex f2bws2 = Resonance::BreitWignerDWave(s2, _f2mass, _f2width,_mpi0,_mpi0);
Complex f2bws3 = Resonance::BreitWignerDWave(s3, _f2mass, _f2width,_mpi0,_mpi0);
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);
+ F1+=Complex(_f2coup*( 0.5*(s3-s2)*f2bws1-Dfact2+Dfact3));
+ F2+=Complex(_f2coup*( 0.5*(s3-s1)*f2bws2-Dfact1+Dfact3));
+ F3+=Complex(_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;
+ F1+=Complex(_f2coup*0.5*(s3-s2)*f2bws1);
+ F2-=Complex(_f2coup*Dfact1);
+ F3-=Complex(_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;
+ F2+=Complex(_f2coup*0.5*(s3-s1)*f2bws2);
+ F1-=Complex(_f2coup*Dfact2);
+ F3+=Complex(_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;
+ F3+=Complex(-_f2coup*0.5*(s1-s2)*f2bws3);
+ F1+=Complex(_f2coup*Dfact3);
+ F2+=Complex(_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;
}
}
// calculate the pi0 pi0 pi+ factor
else if(imode==2) {
// identical particle factors
fact = 1./sqrt(2.);
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3];
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix] = Resonance::BreitWignerPWave(s1,_rhomass[ix], _rhowidth[ix],_mpic,_mpi0);
rhos2bw[ix] = Resonance::BreitWignerPWave(s2,_rhomass[ix], _rhowidth[ix],_mpic,_mpi0);
}
Complex f0bw = Resonance::BreitWignerSWave(s3, _f0mass, _f0width,_mpi0,_mpi0);
Complex sigbw = Resonance::BreitWignerSWave(s3,_sigmamass,_sigmawidth,_mpi0,_mpi0);
Complex f2bw = Resonance::BreitWignerDWave(s3, _f2mass, _f2width,_mpi0,_mpi0);
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]);
+ F1+=Complex(Dfact1*_rhocoupD[ix]*rhos2bw[ix]);
+ F2+=Complex(Dfact2*_rhocoupD[ix]*rhos1bw[ix]);
+ F3+=Complex(_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;
+ F3-=Complex(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];
+ F2+=Complex(Dfact2*_rhocoupD[ires]*rhos1bw[ires]);
+ F3+=Complex(_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];
+ F1+=Complex(Dfact1*_rhocoupD[ires]*rhos2bw[ires]);
+ F3-=Complex(_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;
+ F1+=Dfact3;
+ F2+=Dfact3;
+ F3-=Complex(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(imode==3||imode==4) {
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3];
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix] = Resonance::BreitWignerPWave(s1,_rhomass[ix], _rhowidth[ix],_mpic,_mpi0);
rhos2bw[ix] = Resonance::BreitWignerPWave(s2,_rhomass[ix], _rhowidth[ix],_mpic,_mpi0);
}
Complex f0bw = Resonance::BreitWignerSWave(s3, _f0mass, _f0width,_mpic,_mpic);
Complex sigbw = Resonance::BreitWignerSWave(s3,_sigmamass,_sigmawidth,_mpic,_mpic);
Complex f2bw = Resonance::BreitWignerDWave(s3, _f2mass, _f2width,_mpic,_mpic);
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]);
+ F1+=Complex(Dfact1*_rhocoupD[ix]*rhos2bw[ix]);
+ F2+=Complex(Dfact2*_rhocoupD[ix]*rhos1bw[ix]);
+ F3+=Complex(_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;
+ F3-=Complex(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];
+ F2+=Complex(Dfact2*_rhocoupD[ires]*rhos1bw[ires]);
+ F3+=Complex(_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];
+ F1+=Complex(Dfact1*_rhocoupD[ires]*rhos2bw[ires]);
+ F3-=Complex(_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;
+ F3-=Complex(0.5*_f2coup*(s1-s2)*f2bw);
}
else if(ichan==8) {
F1+=2./3.*_f0coup*f0bw;
F2+=2./3.*_f0coup*f0bw;
}
}
else if(imode==5) {
// identical particle factors
fact = 1./sqrt(2.);
// compute the breit wigners we need
Complex rhos1bw[3],rhos2bw[3];
for(unsigned int ix=0,N=max(_rhocoupP.size(),_rhocoupD.size());ix<N;++ix) {
rhos1bw[ix] = Resonance::BreitWignerPWave(s1,_rhomass[ix], _rhowidth[ix],_mpic,_mpic);
rhos2bw[ix] = Resonance::BreitWignerPWave(s2,_rhomass[ix], _rhowidth[ix],_mpic,_mpic);
}
Complex f0bws1 = Resonance::BreitWignerSWave(s1, _f0mass, _f0width,_mpic,_mpic);
Complex sigbws1 = Resonance::BreitWignerSWave(s1,_sigmamass,_sigmawidth,_mpic,_mpic);
Complex f2bws1 = Resonance::BreitWignerDWave(s1, _f2mass, _f2width,_mpic,_mpic);
Complex f0bws2 = Resonance::BreitWignerSWave(s2, _f0mass, _f0width,_mpic,_mpic);
Complex sigbws2 = Resonance::BreitWignerSWave(s2,_sigmamass,_sigmawidth,_mpic,_mpic);
Complex f2bws2 = Resonance::BreitWignerDWave(s2, _f2mass, _f2width,_mpic,_mpic);
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-=Complex(Dfact1*_rhocoupD[ix]*rhos2bw[ix]);
F2-=Complex(Dfact2*_rhocoupD[ix]*rhos1bw[ix]);
F3-=Complex(_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+=Complex(_f2coup*(0.5*(s3-s2)*f2bws1-sfact2));
F2+=Complex(_f2coup*(0.5*(s3-s1)*f2bws2-sfact1));
F3+=Complex(_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-=Complex(Dfact2*_rhocoupD[ires]*rhos1bw[ires]);
F3-=Complex(_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-=Complex(Dfact1*_rhocoupD[ires]*rhos2bw[ires]);
F3+=Complex(_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+=Complex(_f2coup*0.5*(s3-s2)*f2bws1);
F2-=Complex(_f2coup*sfact1);
F3-=Complex(_f2coup*sfact1);
}
else if(ichan==9) {
complex<Energy2> sfact2 = 1./18.*(4.*_mpic*_mpic-s2)*(q2+s2-_mpic*_mpic)/s2*f2bws2;
F1-=Complex(_f2coup*sfact2);
F2+=Complex(_f2coup*0.5*(s3-s1)*f2bws2);
F3+=Complex(_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;
}
}
else {
throw Exception() << "ThreePionCLEOCurrent Unknown Decay" << imode
<< Exception::abortnow;
}
F1 *= fact;
F2 *= fact;
F3 *= fact;
}
// complete the construction of the decay mode for integration
bool ThreePionCLEOCurrent::createMode(int icharge, tcPDPtr resonance,
FlavourInfo flavour,
unsigned int imode,PhaseSpaceModePtr mode,
unsigned int iloc,int ires,
PhaseSpaceChannel phase, Energy upp ) {
// check the charge and resonance
if(imode<=1||imode==3||imode==4) {
if(icharge!=0) return false;
if(resonance && resonance->id()!=ParticleID::a_10) return false;
}
else if(imode==2||imode==5) {
if(abs(icharge)!=3) return false;
if(resonance && abs(resonance->id())!=ParticleID::a_1plus) return false;
}
// check the total isospin
if(flavour.I!=IsoSpin::IUnknown) {
if(flavour.I!=IsoSpin::IOne) return false;
}
// check I_3
if(flavour.I3!=IsoSpin::I3Unknown) {
switch(flavour.I3) {
case IsoSpin::I3Zero:
if(imode==2||imode==5) return false;
break;
case IsoSpin::I3One:
if((imode!=2&&imode!=5) || icharge ==-3) return false;
break;
case IsoSpin::I3MinusOne:
if((imode!=2&&imode!=5) || icharge == 3) return false;
break;
default:
return false;
}
}
if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return false;
if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return false;
if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return false;
// get the particles and check the masses
int iq(0),ia(0);
tPDVector extpart=particles(1,imode,iq,ia);
Energy min(ZERO);
for(unsigned int ix=0;ix<extpart.size();++ix) min+=extpart[ix]->massMin();
if(min>upp) return false;
_maxmass=max(_maxmass,upp);
// pointers to the particles we need
tPDPtr a1m = getParticleData(ParticleID::a_1minus);
tPDPtr a10 = getParticleData(ParticleID::a_10);
// the different rho resonances
tPDPtr rhom[3] = {getParticleData(-213),getParticleData(-100213),getParticleData(-30213)};
if(icharge==3) {
for(unsigned int ix=0;ix<3;++ix) rhom[ix]=rhom[ix]->CC();
a1m = a1m->CC();
}
tPDPtr rho0[3] = {getParticleData(113),getParticleData(100113),getParticleData(30113)};
// the sigma
tPDPtr sigma = getParticleData(9000221);
// the f_2
tPDPtr f2=getParticleData(225);
// the f_0
tPDPtr f0=getParticleData(10221);
assert(f2 && f0 && sigma);
// a0 -> pi0 pi0 pi0
if(imode<=1) {
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;
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+1,
ires+2,iloc+2,ires+2,iloc+3));
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+2,
ires+2,iloc+1,ires+2,iloc+3));
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+3,
ires+2,iloc+1,ires+2,iloc+2));
}
}
// decay mode a_1- -> pi0 pi0 pi-
else if(imode==2) {
for(unsigned int ix=0;ix<3;++ix) {
// first rho+ channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rhom[ix],ires+1,iloc+1,
ires+2,iloc+2,ires+2,iloc+3));
// second rho+ channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rhom[ix],ires+1,iloc+2,
ires+2,iloc+1,ires+2,iloc+3));
}
// I=0 channels
for(unsigned int iy=0;iy<3;++iy) {
tPDPtr temp;
if(iy==0) temp = sigma;
else if(iy==1) temp = f2;
else if(iy==2) temp = f0;
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+3,
ires+2,iloc+1,ires+2,iloc+2));
}
}
// decay mode a_10 -> pi+ pi- pi0
else if(imode==3||imode==4) {
// rho modes
for(unsigned int ix=0;ix<3;++ix) {
// first rho channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,rhom[ix],ires+1,iloc+1,
ires+2,iloc+2,ires+2,iloc+3));
// second channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,rhom[ix],ires+1,iloc+2,
ires+2,iloc+1,ires+2,iloc+3));
}
// I=0 channels
for(unsigned int iy=0;iy<3;++iy) {
tPDPtr temp;
if(iy==0) temp = sigma;
else if(iy==1) temp = f2;
else if(iy==2) temp = f0;
mode->addChannel((PhaseSpaceChannel(phase),ires,a10,ires+1,temp,ires+1,iloc+3,
ires+2,iloc+1,ires+2,iloc+2));
}
}
else if(imode==5) {
for(unsigned int ix=0;ix<3;++ix) {
// the neutral rho channels
// first channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rho0[ix],ires+1,iloc+1,
ires+2,iloc+2,ires+2,iloc+3));
// interchanged channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,rho0[ix],ires+1,iloc+2,
ires+2,iloc+1,ires+2,iloc+3));
}
for(unsigned int iy=0;iy<3;++iy) {
tPDPtr temp;
if(iy==0) temp = sigma;
else if(iy==1) temp = f2;
else if(iy==2) temp = f0;
// first channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+1,
ires+2,iloc+2,ires+2,iloc+3));
// interchanged channel
mode->addChannel((PhaseSpaceChannel(phase),ires,a1m,ires+1,temp,ires+1,iloc+2,
ires+2,iloc+1,ires+2,iloc+3));
}
}
// reset the integration parameters
for(unsigned int iy=0;iy<_rhomass.size();++iy) {
mode->resetIntermediate(rho0[iy],_rhomass[iy],_rhowidth[iy]);
mode->resetIntermediate(rhom[iy],_rhomass[iy],_rhowidth[iy]);
}
mode->resetIntermediate(sigma,_sigmamass,_sigmawidth);
mode->resetIntermediate(f2,_f2mass,_f2width);
mode->resetIntermediate(f0,_f0mass,_f0width);
mode->resetIntermediate(a10,_a1mass,_a1width);
mode->resetIntermediate(a10,_a1mass,_a1width);
return true;
}
void ThreePionCLEOCurrent::dataBaseOutput(ofstream & output,bool header,
bool create) const {
if(header){output << "update decayers set parameters=\"";}
if(create) {
output << "create Herwig::ThreePionCLEOCurrent " << name()
<< " HwWeakCurrents.so\n";
}
for(unsigned int ix=0;ix<_rhomass.size();++ix) {
if(ix<2) {
output << "newdef " << name() << ":RhoMasses " << ix
<< " " << _rhomass[ix]/MeV << "\n";
}
else {
output << "insert " << name() << ":RhoMasses " << ix
<< " " << _rhomass[ix]/MeV << "\n";
}
}
for(unsigned int ix=0;ix<_rhowidth.size();++ix) {
if(ix<2) {
output << "newdef " << name() << ":RhoWidths " << ix
<< " " << _rhowidth[ix]/MeV << "\n";
}
else {
output << "insert " << name() << ":RhoWidths " << ix
<< " " << _rhowidth[ix]/MeV << "\n";
}
}
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";
output << "newdef " << name() << ":a1Mass " << _a1mass/GeV << "\n";
output << "newdef " << name() << ":a1Width " <<_a1width /GeV << "\n";
output << "newdef " << name() << ":KaonMass " << _mK/GeV << "\n";
output << "newdef " << name() << ":KStarMass " << _mKstar/GeV << "\n";
output << "newdef " << name() << ":KaonCoupling " << _gammk << "\n";
output << "newdef " << name() << ":Fpi " << _fpi/MeV << "\n";
output << "newdef " << name() << ":a1WidthOption " << _a1opt << "\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";
}
}
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() << ":Initializea1 " <<_initializea1 << "\n";
for(unsigned int ix=0;ix<_a1runwidth.size();++ix) {
if(ix<200) {
output << "newdef " << name() << ":a1RunningWidth " << ix
<< " " << _a1runwidth[ix]/MeV << "\n";
}
else {
output << "insert " << name() << ":a1RunningWidth " << ix
<< " " << _a1runwidth[ix]/MeV << "\n";
}
}
for(unsigned int ix=0;ix<_a1runq2.size();++ix) {
if(ix<200) {
output << "newdef " << name() << ":a1RunningQ2 " << ix
<< " " << _a1runq2[ix]/MeV2 << "\n";
}
else {
output << "insert " << name() << ":a1RunningQ2 " << ix
<< " " << _a1runq2[ix]/MeV2 << "\n";
}
}
WeakCurrent::dataBaseOutput(output,false,false);
if(header) output << "\n\" where BINARY ThePEGName=\""
<< fullName() << "\";" << endl;
}
void ThreePionCLEOCurrent::doinitrun() {
// set up the running a_1 width
inita1Width(0);
WeakCurrent::doinitrun();
}
void ThreePionCLEOCurrent::doupdate() {
WeakCurrent::doupdate();
// update running width if needed
if ( !touched() ) return;
if(_maxmass!=_maxcalc) inita1Width(-1);
}
Energy ThreePionCLEOCurrent::a1width(Energy2 q2) const {
Energy output;
if(_a1opt) output=(*_a1runinter)(q2);
else {
double gam(0.);
if(q2<0.1753*GeV2) {
gam =0.;
}
else if(q2<0.823*GeV2) {
double p=q2/GeV2-0.1753;
gam = 5.80900*p*sqr(p)*(1.-3.00980*p+4.57920*sqr(p));
}
else {
double p=q2/GeV2;
gam = -13.91400+27.67900*p-13.39300*sqr(p)
+3.19240*sqr(p)*p-0.10487*sqr(sqr(p));
}
if(q2<0.1676*GeV2) {
gam+=0.;
}
else if(q2<0.823*GeV2) {
double p=q2/GeV2-0.1676;
gam+= 6.28450*p*sqr(p)*(1.-2.95950*p+4.33550*sqr(p));
}
else {
double p=q2/GeV2;
gam+= -15.41100+32.08800*p-17.66600*sqr(p)
+4.93550*sqr(p)*p-0.37498*sqr(sqr(p));
}
Energy mkst=0.894*GeV,mk=0.496*GeV;
Energy2 mk1sq=sqr(mkst+mk), mk2sq=sqr(mkst-mk);
double c3pi=sqr(0.2384),ckst=sqr(4.7621)*c3pi;
gam*=c3pi;
if(q2>mk1sq) gam+=0.5*ckst*sqrt((q2-mk1sq)*(q2-mk2sq))/q2;
gam = gam*_a1width*_a1mass/GeV2/1.331/0.814/1.0252088;
output = gam*GeV2/sqrt(q2);
}
return output;
}
double
ThreePionCLEOCurrent::threeBodyMatrixElement(const int iopt, const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy,
const Energy, const Energy) const {
Energy p1[5],p2[5],p3[5];
Energy2 p1sq, p2sq, p3sq;
Energy q=sqrt(q2);
Energy2 mpi2c=_mpic*_mpic;
Energy2 mpi20=_mpi0*_mpi0;
// construct the momenta for the 2 neutral 1 charged mode
Complex F1,F2,F3;
if(iopt==0) {
// construct the momenta of the decay products
p1[0] = 0.5*(q2+mpi20-s1)/q; p1sq=p1[0]*p1[0]; p1[4]=sqrt(p1sq-mpi20);
p2[0] = 0.5*(q2+mpi20-s2)/q; p2sq=p2[0]*p2[0]; p2[4]=sqrt(p2sq-mpi20);
p3[0] = 0.5*(q2+mpi2c-s3)/q; p3sq=p3[0]*p3[0]; p3[4]=sqrt(p3sq-mpi2c);
// take momentum of 1 parallel to z axis
p1[1]=ZERO;p1[2]=ZERO;p1[3]=p1[4];
// construct 2
double cos2 = 0.5*(p1sq+p2sq-p3sq-2.*mpi20+mpi2c)/p1[4]/p2[4];
p2[1] = p2[4]*sqrt(1.-cos2*cos2); p2[2]=ZERO; p2[3]=-p2[4]*cos2;
// construct 3
double cos3 = 0.5*(p1sq-p2sq+p3sq-mpi2c)/p1[4]/p3[4];
p3[1] =-p3[4]*sqrt(1.-cos3*cos3); p3[2]=ZERO; p3[3]=-p3[4]*cos3;
// calculate the form factors
CLEOFormFactor(1,-1,q2,s1,s2,s3,F1,F2,F3);
}
// construct the momenta for the 3 charged mode
else {
// construct the momenta of the decay products
p1[0] = 0.5*(q2+mpi2c-s1)/q; p1sq=p1[0]*p1[0]; p1[4]=sqrt(p1sq-mpi2c);
p2[0] = 0.5*(q2+mpi2c-s2)/q; p2sq=p2[0]*p2[0]; p2[4]=sqrt(p2sq-mpi2c);
p3[0] = 0.5*(q2+mpi2c-s3)/q; p3sq=p3[0]*p3[0]; p3[4]=sqrt(p3sq-mpi2c);
// take momentum of 1 parallel to z axis
p1[1]=ZERO;p1[2]=ZERO;p1[3]=p1[4];
// construct 2
double cos2 = 0.5*(p1sq+p2sq-p3sq-mpi2c)/p1[4]/p2[4];
p2[1] = p2[4]*sqrt(1.-cos2*cos2); p2[2]=ZERO; p2[3]=-p2[4]*cos2;
// construct 3
double cos3 = 0.5*(p1sq-p2sq+p3sq-mpi2c)/p1[4]/p3[4];
p3[1] =-p3[4]*sqrt(1.-cos3*cos3); p3[2]=ZERO; p3[3]=-p3[4]*cos3;
// calculate the form factors
CLEOFormFactor(0,-1,q2,s1,s2,s3,F1,F2,F3);
}
// construct a vector with the current
complex<Energy> current[4];
for(unsigned int ix=0;ix<4;++ix)
current[ix] = F1*(p2[ix]-p3[ix])-F2*(p3[ix]-p1[ix])+F3*(p1[ix]-p2[ix]);
complex<Energy2> dot1=current[0]*conj(current[0]);
for(unsigned int ix=1;ix<4;++ix) dot1-=current[ix]*conj(current[ix]);
complex<Energy2> dot2=current[0]*q;
return(-dot1+dot2*conj(dot2)/q2).real() / sqr(_rhomass[0]);
}
// the hadronic currents
vector<LorentzPolarizationVectorE>
ThreePionCLEOCurrent::current(tcPDPtr resonance,
FlavourInfo flavour,
const int imode, const int ichan, Energy & scale,
const tPDVector & ,
const vector<Lorentz5Momentum> & momenta,
DecayIntegrator::MEOption) const {
useMe();
// check the isospin
if(flavour.I!=IsoSpin::IUnknown && flavour.I!=IsoSpin::IOne)
return vector<LorentzPolarizationVectorE>();
// check I_3
if(flavour.I3!=IsoSpin::I3Unknown) {
switch(flavour.I3) {
case IsoSpin::I3Zero:
if(imode==2||imode==5) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3One:
if(imode!=2&&imode!=5) return vector<LorentzPolarizationVectorE>();
break;
case IsoSpin::I3MinusOne:
if(imode!=2&&imode!=5) return vector<LorentzPolarizationVectorE>();
break;
default:
return vector<LorentzPolarizationVectorE>();
}
}
if(flavour.strange != Strangeness::Unknown and flavour.strange != Strangeness::Zero) return vector<LorentzPolarizationVectorE>();
if(flavour.charm != Charm::Unknown and flavour.charm != Charm::Zero ) return vector<LorentzPolarizationVectorE>();
if(flavour.bottom != Beauty::Unknown and flavour.bottom !=Beauty::Zero ) return vector<LorentzPolarizationVectorE>();
// calculate q2,s1,s2,s3
Lorentz5Momentum q;
for(unsigned int ix=0;ix<momenta.size();++ix)
q+=momenta[ix];
q.rescaleMass();
scale=q.mass();
Energy2 q2=q.mass2();
Energy2 s1 = (momenta[1]+momenta[2]).m2();
Energy2 s2 = (momenta[0]+momenta[2]).m2();
Energy2 s3 = (momenta[0]+momenta[1]).m2();
// form factors
Complex F1(0.), F2(0.), F3(0.);
CLEOFormFactor(imode,ichan,q2,s1,s2,s3,F1,F2,F3);
// change sign of the F2 term
F2 =- F2;
// prefactor
complex<InvEnergy> a1fact = _fact;
if(!resonance) a1fact *= a1BreitWigner(q2);
// current
LorentzPolarizationVectorE vect = q.mass()*a1fact*
((F2-F1)*momenta[2] + (F1-F3)*momenta[1] + (F3-F2)*momenta[0]);
// scalar piece
Complex dot=(vect*q)/q2;
vect -= dot*q;
// return the answer
return vector<LorentzPolarizationVectorE>(1,vect);
}
bool ThreePionCLEOCurrent::accept(vector<int> id) {
if(id.size()!=3) return false;
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID::piplus) continue;
else if(id[ix]==ParticleID::piminus) continue;
else if(id[ix]==ParticleID::pi0) continue;
return false;
}
return true;
}
unsigned int ThreePionCLEOCurrent::decayMode(vector<int> id) {
if(id.size()!=3) return -1;
int npip(0),npim(0),npi0(0);
for(unsigned int ix=0;ix<id.size();++ix) {
if(id[ix]==ParticleID::piplus) ++npip;
else if(id[ix]==ParticleID::piminus) ++npim;
else if(id[ix]==ParticleID::pi0) ++npi0;
}
if (npi0==3) return 0;
else if( (npip==1&&npi0==2) || (npim==1&&npi0==2) ) return 2;
else if( npi0==1 && npip==1 && npim==1 ) return 3;
else if( (npip==2&&npim==1) || (npim==2&&npip==1) ) return 5;
else return -1;
}
tPDVector ThreePionCLEOCurrent::particles(int icharge, unsigned int imode,int,int) {
tPDVector extpart(3);
if(imode==0||imode==1) {
extpart[0]=getParticleData(ParticleID::pi0);
extpart[1]=getParticleData(ParticleID::pi0);
extpart[2]=getParticleData(ParticleID::pi0);
}
else if(imode==2) {
extpart[0]=getParticleData(ParticleID::pi0);
extpart[1]=getParticleData(ParticleID::pi0);
extpart[2]=getParticleData(ParticleID::piminus);
}
else if(imode==3||imode==4) {
extpart[0]=getParticleData(ParticleID::piplus);
extpart[1]=getParticleData(ParticleID::piminus);
extpart[2]=getParticleData(ParticleID::pi0);
}
else if(imode==5) {
extpart[0]=getParticleData(ParticleID::piminus);
extpart[1]=getParticleData(ParticleID::piminus);
extpart[2]=getParticleData(ParticleID::piplus);
}
else
assert(false);
// conjugate the particles if needed
if(icharge==3) {
for(unsigned int ix=0;ix<3;++ix) {
if(extpart[ix]->CC()) extpart[ix]=extpart[ix]->CC();
}
}
// return the answer
return extpart;
}
diff --git a/Hadronization/SpinHadronizer.h b/Hadronization/SpinHadronizer.h
--- a/Hadronization/SpinHadronizer.h
+++ b/Hadronization/SpinHadronizer.h
@@ -1,191 +1,191 @@
// -*- C++ -*-
#ifndef Herwig_SpinHadronizer_H
#define Herwig_SpinHadronizer_H
//
// This is the declaration of the SpinHadronizer class.
//
#include "ThePEG/Handlers/StepHandler.h"
namespace Herwig {
using namespace ThePEG;
/**
* The SpinHadronizer class is designed to be used as a post-hadronization handler to
* give a simple model of spin transfer between the perturbative and non-perturbative
* stages.
*
* @see \ref SpinHadronizerInterfaces "The interfaces"
* defined for SpinHadronizer.
*/
class SpinHadronizer: public StepHandler {
public:
/**
* The default constructor.
*/
SpinHadronizer() : omegaHalf_(2./3.), omegaThreeHalf_(0.2),
minFlav_(3), maxFlav_(5), debug_(false), qPol_(6,make_pair(0.,0.))
{}
public:
/** @name Virtual functions required by the StepHandler class. */
//@{
/**
* The main function called by the EventHandler class to
* perform a step. Given the current state of an Event, this function
* performs the event generation step and includes the result in a new
* Step object int the Event record.
* @param eh the EventHandler in charge of the Event generation.
* @param tagged if not empty these are the only particles which should
* be considered by the StepHandler.
* @param hint a Hint object with possible information from previously
* performed steps.
* @throws Veto if the StepHandler requires the current step to be discarded.
* @throws Stop if the generation of the current Event should be stopped
* after this call.
* @throws Exception if something goes wrong.
*/
virtual void handle(EventHandler & eh, const tPVector & tagged,
const Hint & hint);
//@}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Functions to calculate the spins
*/
//@{
/**
* Calculate the spin of a baryon
*/
void baryonSpin(tPPtr baryon);
/**
* Calculate the spin of a meson
*/
void mesonSpin(tPPtr meson);
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- SpinHadronizer & operator=(const SpinHadronizer &);
+ SpinHadronizer & operator=(const SpinHadronizer &) = delete;
private:
/**
* Parameters
*/
//@{
/**
* Falk-Peskin \f$\omega_\frac12\f$ parameter
*/
double omegaHalf_;
/**
* Falk-Peskin \f$\omega_\frac32\f$ parameter
*/
double omegaThreeHalf_;
/**
* Minimum quark flavour
*/
unsigned int minFlav_;
/**
* Maximum quark flavour
*/
unsigned int maxFlav_;
/**
* Print out debugging info
*/
bool debug_;
/**
* Polarization of the quarks
*/
vector<pair<double,double> > qPol_;
//@}
};
}
#endif /* Herwig_SpinHadronizer_H */
diff --git a/MatrixElement/Matchbox/Utility/DensityOperator.h b/MatrixElement/Matchbox/Utility/DensityOperator.h
--- a/MatrixElement/Matchbox/Utility/DensityOperator.h
+++ b/MatrixElement/Matchbox/Utility/DensityOperator.h
@@ -1,248 +1,248 @@
// -*- C++ -*-
//
// DensityOperator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_DensityOperator_H
#define Herwig_DensityOperator_H
//
// This is the declaration of the DensityOperator class.
//
#include "ThePEG/Handlers/HandlerBase.h"
#include "Herwig/MatrixElement/Matchbox/Utility/ColourBasis.h"
#include <tuple>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
namespace Herwig {
using namespace ThePEG;
typedef boost::numeric::ublas::vector<Complex> CVector;
/**
* Here is the documentation of the DensityOperator class.
*
* @see \ref DensityOperatorInterfaces "The interfaces"
* defined for DensityOperator.
*/
class DensityOperator: public HandlerBase {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DensityOperator();
/**
* The destructor.
*/
virtual ~DensityOperator();
//@}
public:
/**
* Clears theDensityOperatorMap.
*/
void clear();
/**
* Prepare for the given sub process.
*/
void prepare(const cPDVector&);
/**
* Fill the density operator for the given hard subprocess, summing over all
* helicity configurations.
*/
void fill(const Ptr<MatchboxXComb>::ptr,
const cPDVector&, const vector<Lorentz5Momentum>& momenta);
/**
* Evolve the density operator, by
* M_{n+1} = -\sum_{i,k}{-4*pi*alpha_s/Ti2*V_{ij,k} T_{i,n}M_nT_{k,n}^\dag},
* see arXiv:1206.0180 eq. (5), note that the pi*pj factor is assumed to be
* included in V_{ij,k}.
*/
void evolve(const map<pair<size_t,size_t>,Complex>& Vijk,
const cPDVector& before,
const cPDVector& after,
const map<std::tuple<size_t,size_t,size_t>,map<size_t,size_t> >& emissionsMap,
const bool splitAGluon,
const bool initialGluonSplitting);
/**
* Calculate the colour matrix element correction.
* -(1+delta(i,gluon))/Ti^2 Tr(Sn+1 Ti Mn Tk^dagger)/Tr(Sn Mn)
* where the bracket in front compensates for the gluon symmetry factor,
* Ti^2 is C_f or C_a, Sn+1 is the matrix of scalar products, and
* Ti is the radiation matrix.
* The first arg contains (emitter index, spectator index, emission pid)
*
*/
double colourMatrixElementCorrection(const std::tuple<size_t,size_t,long>& ikemission,
const cPDVector& particles);
/**
* Checking colour conservation for the colour matrix element corrections.
*/
void colourConservation(const cPDVector& particles);
/**
* Get the colour basis.
*/
Ptr<ColourBasis>::tptr colourBasis() { return theColourBasis; }
/**
* Get the colour basis.
*/
const Ptr<ColourBasis>::tptr colourBasis() const { return theColourBasis; }
/**
* Set the colour basis.
*/
void colourBasis(Ptr<ColourBasis>::ptr ptr) { theColourBasis = ptr; }
/**
* Get the correlator map.
*/
const map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double>& correlatorMap() const {
return theCorrelatorMap;
}
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Number of colours used in colourNorm.
*/
double Nc;
/**
* QCD vertex normalization.
*/
double TR;
/**
* Normalization of colour charges \mathbf{T}_{ij}^2.
*/
double colourNorm(const cPDPtr particle);
/**
* Fast evaluation of Tij*Mn, where a Tij is the matrix from ColourBasis::charge,
* which is a sparse matrix, and Mn is the density operator, a dense matrix.
*
*/
matrix<Complex> prodSparseDense(const compressed_matrix<double>&,
const matrix<Complex>&);
/**
* Fast evaluation of TijMn*Tkdagger, where a TijMn is the result from the method
* prodSparseDense, a dense matrix, and Tkdagger is the transponse conjugate of
* the matrix from ColourBasis::charge, a sparse matrix.
*
*/
matrix<Complex> prodDenseSparse(const matrix<Complex>&,
const compressed_matrix<double>&);
/**
* Boosts a vector of momenta to the rest frame of the initial pair
* of particles (the first 2 elements of the argument vector). Returns
* the boosted vectors
*/
vector<Lorentz5Momentum> boostToRestFrame(const vector<Lorentz5Momentum>& momenta);
/**
* Boosts a vector of momenta to the rest frame of the initial pair
*/
bool compareMomentum(const Lorentz5Momentum& p, const Lorentz5Momentum& q);
/**
* Mapping of colour structures to density operator matrices.
*
*/
map<vector<PDT::Colour>,matrix<Complex> > theDensityOperatorMap;
/**
* Mapping of colour structures and legs to colour correlators.
*/
map<pair<vector<PDT::Colour>,pair<size_t,size_t> >,double> theCorrelatorMap;
/**
* A map from the hard subprocess particles to a map of amplitude colour
* basis order to the normal ordered colour basis.
*/
map<cPDVector, map<size_t,size_t> > theColourBasisToColourBasisMap;
/**
* Colour basis used.
*/
Ptr<ColourBasis>::ptr theColourBasis;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- DensityOperator & operator=(const DensityOperator &);
+ DensityOperator & operator=(const DensityOperator &) = delete;
};
}
#endif /* Herwig_DensityOperator_H */
diff --git a/Models/DarkMatter/DMDMMediatorVertex.h b/Models/DarkMatter/DMDMMediatorVertex.h
--- a/Models/DarkMatter/DMDMMediatorVertex.h
+++ b/Models/DarkMatter/DMDMMediatorVertex.h
@@ -1,113 +1,113 @@
// -*- C++ -*-
#ifndef Herwig_DMDMMediatorVertex_H
#define Herwig_DMDMMediatorVertex_H
//
// This is the declaration of the DMDMMediatorVertex class.
//
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* The DMDMMediatorVertex class implements the coupling of the dark matter to the mediator
*
* @see \ref DMDMMediatorVertexInterfaces "The interfaces"
* defined for DMDMMediatorVertex.
*/
class DMDMMediatorVertex: public FFVVertex {
public:
/**
* The default constructor.
*/
DMDMMediatorVertex();
/**
* Calculate the couplings.
* @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
* @param part1 The ParticleData pointer for the first particle.
* @param part2 The ParticleData pointer for the second particle.
* @param part3 The ParticleData pointer for the third particle.
*/
virtual void setCoupling(Energy2 q2,tcPDPtr part1,tcPDPtr part2,tcPDPtr part3);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- DMDMMediatorVertex & operator=(const DMDMMediatorVertex &);
+ DMDMMediatorVertex & operator=(const DMDMMediatorVertex &) = delete;
private:
/**
* DM coupling to the dark mediator
*/
double cDMmed_;
};
}
#endif /* Herwig_DMDMMediatorVertex_H */
diff --git a/Models/DarkMatter/MEDM2Mesons.cc b/Models/DarkMatter/MEDM2Mesons.cc
--- a/Models/DarkMatter/MEDM2Mesons.cc
+++ b/Models/DarkMatter/MEDM2Mesons.cc
@@ -1,266 +1,266 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the MEDM2Mesons class.
//
#include "MEDM2Mesons.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "Herwig/Models/General/BSMModel.h"
using namespace Herwig;
typedef LorentzVector<complex<InvEnergy> > LorentzPolarizationVectorInvE;
MEDM2Mesons::MEDM2Mesons() : cDMmed_(0.), cSMmed_({0.0,0.0,0.0})
{}
Energy2 MEDM2Mesons::scale() const {
return sHat();
}
unsigned int MEDM2Mesons::orderInAlphaS() const {
return 0;
}
unsigned int MEDM2Mesons::orderInAlphaEW() const {
return 0;
}
IBPtr MEDM2Mesons::clone() const {
return new_ptr(*this);
}
IBPtr MEDM2Mesons::fullclone() const {
return new_ptr(*this);
}
void MEDM2Mesons::doinit() {
// make sure the current got initialised
current_->init();
// max energy
Energy Emax = generator()->maximumCMEnergy();
// loop over the modes
int nmode=0;
for(unsigned int imode=0;imode<current_->numberOfModes();++imode) {
// get the external particles for this mode
int iq(0),ia(0);
tPDVector out = current_->particles(0,imode,iq,ia);
current_->decayModeInfo(imode,iq,ia);
if(iq==2&&ia==-2) continue;
PhaseSpaceModePtr mode = new_ptr(PhaseSpaceMode(incomingA_,out,1.,
incomingB_,Emax));
PhaseSpaceChannel channel(mode);
if(!current_->createMode(0,tcPDPtr(), FlavourInfo(), imode,mode,0,-1,channel,Emax)) continue;
modeMap_[imode] = nmode;
addMode(mode);
++nmode;
}
// cast the model
Ptr<BSMModel>::ptr model = dynamic_ptr_cast<Ptr<BSMModel>::ptr>(generator()->standardModel());
bool foundDM(false),foundU(false),foundD(false),foundS(false);
// find the vertices we need and extract the couplings
for(unsigned int ix = 0; ix < model->numberOfVertices(); ++ix ) {
VertexBasePtr vertex = model->vertex(ix);
if(vertex->getNpoint()!=3) continue;
for(unsigned int iloc = 0;iloc < 3; ++iloc) {
vector<long> ext = vertex->search(iloc, Mediator_->id());
if(ext.empty()) continue;
for(unsigned int ioff=0;ioff<ext.size();ioff+=3) {
if(iloc!=2) assert(false);
if((ext[ioff]==incomingA_->id() && ext[ioff+1]==incomingB_->id()) ||
(ext[ioff]==incomingB_->id() && ext[ioff+1]==incomingA_->id())) {
foundDM = true;
vertex->setCoupling(sqr(Emax),incomingA_,incomingB_,Mediator_);
cDMmed_ = vertex->norm();
}
else if(abs(ext[ioff])==1 && abs(ext[ioff+1])==1 && ext[ioff]==-ext[ioff+1]) {
foundD = true;
vertex->setCoupling(sqr(Emax),getParticleData(1),getParticleData(-1),Mediator_);
cSMmed_[0] = vertex->norm();
}
else if(abs(ext[ioff])==2 && abs(ext[ioff+1])==2 && ext[ioff]==-ext[ioff+1]) {
foundU = true;
vertex->setCoupling(sqr(Emax),getParticleData(2),getParticleData(-2),Mediator_);
cSMmed_[1] = vertex->norm();
}
else if(abs(ext[ioff])==3 && abs(ext[ioff+1])==3 && ext[ioff]==-ext[ioff+1]) {
foundS = true;
vertex->setCoupling(sqr(Emax),getParticleData(3),getParticleData(-3),Mediator_);
cSMmed_[2] = vertex->norm();
}
}
}
}
if(!foundDM) {
throw InitException() << "Cannot find DM coupling in MEDM2Mesons::doinit()";
}
if(!foundD) {
throw InitException() << "Cannot find down quark coupling in MEDM2Mesons::doinit()";
}
if(!foundU) {
throw InitException() << "Cannot find up quark coupling in MEDM2Mesons::doinit()";
}
if(!foundS) {
throw InitException() << "Cannot find strange quark coupling in MEDM2Mesons::doinit()";
}
MEMultiChannel::doinit();
}
void MEDM2Mesons::persistentOutput(PersistentOStream & os) const {
os << current_ << modeMap_ << incomingA_ << incomingB_ << Mediator_ << cDMmed_ << cSMmed_;
}
void MEDM2Mesons::persistentInput(PersistentIStream & is, int) {
is >> current_ >> modeMap_ >> incomingA_ >> incomingB_ >> Mediator_ >> cDMmed_ >> cSMmed_;
}
//The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<MEDM2Mesons,MEMultiChannel>
describeHerwigMEDM2Mesons("Herwig::MEDM2Mesons", "Herwig.so");
void MEDM2Mesons::Init() {
static ClassDocumentation<MEDM2Mesons> documentation
("The MEDM2Mesons class simulates the annhilation of"
" DM particles to mesons at low energy");
static Reference<MEDM2Mesons,WeakCurrent> interfaceWeakCurrent
("WeakCurrent",
"The reference for the decay current to be used.",
&MEDM2Mesons::current_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceIncomingA
("IncomingA",
"First incoming particle",
&MEDM2Mesons::incomingA_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceIncomingB
("IncomingB",
"Second incoming particle",
&MEDM2Mesons::incomingB_, false, false, true, false, false);
static Reference<MEDM2Mesons,ParticleData> interfaceMediator_
("Mediator",
"DM mediator",
&MEDM2Mesons::Mediator_, false, false, true, false, false);
}
double MEDM2Mesons::me2(const int ichan) const {
// compute the incoming current
LorentzPolarizationVectorInvE lepton[2][2];
if(incomingA_->iSpin()==PDT::Spin1Half && incomingB_->iSpin()==PDT::Spin1Half) {
SpinorWaveFunction em_in( meMomenta()[0],mePartonData()[0],incoming);
SpinorBarWaveFunction ep_in( meMomenta()[1],mePartonData()[1],incoming);
vector<SpinorWaveFunction> f1;
vector<SpinorBarWaveFunction> a1;
for(unsigned int ix=0;ix<2;++ix) {
em_in.reset(ix);
f1.push_back(em_in);
ep_in.reset(ix);
a1.push_back(ep_in);
}
// this should be coupling of DM to mediator/ mediator propagator
complex<Energy> mmed = Mediator_->mass();
complex<Energy2> mmed2 = sqr(mmed);
complex<Energy> mwid = Mediator_->width();
complex<Energy2> prop = sHat()-mmed2+Complex(0.,1.)*mmed*mwid;
complex<InvEnergy2> pre = cDMmed_/prop;
for(unsigned ix=0;ix<2;++ix) {
for(unsigned iy=0;iy<2;++iy) {
lepton[ix][iy]= pre*f1[ix].dimensionedWave().vectorCurrent(a1[iy].dimensionedWave());
}
}
}
// TODO think about other spins for the DM
else
assert(false);
// work out the mapping for the hadron vector
int nOut = int(meMomenta().size())-2;
vector<unsigned int> constants(nOut+1);
vector<PDT::Spin > iSpin(nOut);
vector<int> hadrons(nOut);
int itemp(1);
int ix(nOut);
do {
--ix;
iSpin[ix] = mePartonData()[ix+2]->iSpin();
itemp *= iSpin[ix];
constants[ix] = itemp;
hadrons[ix] = mePartonData()[ix+2]->id();
}
while(ix>0);
constants[nOut] = 1;
// calculate the matrix element
me_.reset(ProductionMatrixElement(PDT::Spin1Half,PDT::Spin1Half,iSpin));
// calculate the hadron current
unsigned int imode = current_->decayMode(hadrons);
Energy q = sqrt(sHat());
vector<Lorentz5Momentum> momenta(meMomenta() .begin()+2, meMomenta().end());
tPDVector out = mode(modeMap_.at(imode))->outgoing();
if(ichan<0) iMode(modeMap_.at(imode));
// get the hadronic currents for the I=1 and I=0 components
vector<LorentzPolarizationVectorE>
hadronI0(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::Zero),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronI1(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IOne, IsoSpin::I3Zero,Strangeness::Zero),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
vector<LorentzPolarizationVectorE>
hadronssbar(current_->current(tcPDPtr(), FlavourInfo(IsoSpin::IZero, IsoSpin::I3Zero,Strangeness::ssbar),
imode,ichan,q,out,momenta,DecayIntegrator::Calculate));
// compute the matrix element
vector<unsigned int> ihel(meMomenta().size());
double output(0.);
unsigned int hI0_size = hadronI0.size();
unsigned int hI1_size = hadronI1.size();
unsigned int hss_size = hadronssbar.size();
unsigned int maxsize = max(max(hadronI0.size(),hadronI1.size()),hss_size);
for(unsigned int hhel=0;hhel<maxsize;++hhel) {
// map the index for the hadrons to a helicity state
for(int ix=nOut;ix>0;--ix) {
ihel[ix+1]=(hhel%constants[ix-1])/constants[ix];
}
// loop over the helicities of the incoming particles
for(ihel[1]=0;ihel[1]<2;++ihel[1]){
for(ihel[0]=0;ihel[0]<2;++ihel[0]) {
Complex amp = 0.;
// work on coefficients for the I1 and I0 bits
if(hI0_size != 0 )
- amp += (cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel]));
+ amp += Complex((cSMmed_[0]+cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI0[hhel])));
if(hI1_size !=0)
- amp += (cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel]));
+ amp += Complex((cSMmed_[0]-cSMmed_[1])/sqrt(2.)*(lepton[ihel[0]][ihel[1]].dot(hadronI1[hhel])));
if(hss_size !=0)
- amp += cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel]));
+ amp += Complex(cSMmed_[2]* (lepton[ihel[0]][ihel[1]].dot(hadronssbar[hhel])));
me_(ihel)= amp;
output += std::norm(amp);
}
}
}
// symmetry factors
map<long,int> ncount;
double symmetry(1.);
for(tPDPtr o : out) ncount[o->id()]+=1;
for(map<long,int>::const_iterator it=ncount.begin();it!=ncount.end();++it) {
symmetry *= it->second;
}
// prefactors
output *= 0.25*sqr(pow(sqrt(sHat())/q,int(momenta.size()-2)));
return output/symmetry;
}
void MEDM2Mesons::constructVertex(tSubProPtr) {
}
diff --git a/Models/DarkMatter/MEDM2Mesons.h b/Models/DarkMatter/MEDM2Mesons.h
--- a/Models/DarkMatter/MEDM2Mesons.h
+++ b/Models/DarkMatter/MEDM2Mesons.h
@@ -1,177 +1,177 @@
// -*- C++ -*-
#ifndef Herwig_MEDM2Mesons_H
#define Herwig_MEDM2Mesons_H
//
// This is the declaration of the MEDM2Mesons class.
//
#include "Herwig/MatrixElement/MEMultiChannel.h"
#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
#include "Herwig/MatrixElement/ProductionMatrixElement.h"
namespace Herwig {
using namespace ThePEG;
/**
* The MEDM2Mesons class implements dark matter annhilation via a vector current to
* mesons at low energies
*
* @see \ref MEDM2MesonsInterfaces "The interfaces"
* defined for MEDM2Mesons.
*/
class MEDM2Mesons: public MEMultiChannel {
public:
/**
* The default constructor.
*/
MEDM2Mesons();
public:
/** @name Virtual functions required by the MEBase class. */
//@{
/**
* Return the order in \f$\alpha_S\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaS() const;
/**
* Return the order in \f$\alpha_{EW}\f$ in which this matrix
* element is given.
*/
virtual unsigned int orderInAlphaEW() const;
/**
* Return the scale associated with the last set phase space point.
*/
virtual Energy2 scale() const;
//@}
/**
* Construct the vertex of spin correlations.
*/
virtual void constructVertex(tSubProPtr);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
*/
virtual double me2(const int ichan) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
private:
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- MEDM2Mesons & operator=(const MEDM2Mesons &);
+ MEDM2Mesons & operator=(const MEDM2Mesons &) = delete;
private :
/**
* Hadronic current etc
*/
//@{
/**
* the hadronic current
*/
WeakCurrentPtr current_;
/**
* The matrix element
*/
mutable ProductionMatrixElement me_;
/**
* Map for the modes
*/
map<int,int> modeMap_;
//@}
/**
* DM
*/
//@{
/**
* Incoming Particles
*/
PDPtr incomingA_, incomingB_;
//@}
/**
* DM coupling to the dark mediator
*/
Complex cDMmed_;
/**
* SM couplings to the dark mediator
*/
vector<Complex> cSMmed_;
/**
* DM vector mediator
*/
PDPtr Mediator_;
};
}
#endif /* Herwig_MEDM2Mesons_H */
diff --git a/PDF/HwRemDecayer.h b/PDF/HwRemDecayer.h
--- a/PDF/HwRemDecayer.h
+++ b/PDF/HwRemDecayer.h
@@ -1,747 +1,747 @@
// -*- C++ -*-
//
// HwRemDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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_HwRemDecayer_H
#define HERWIG_HwRemDecayer_H
//
// This is the declaration of the HwRemDecayer class.
//
#include "ThePEG/PDT/RemnantDecayer.h"
#include "ThePEG/Handlers/EventHandler.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/EventRecord/SubProcess.h"
#include "ThePEG/PDF/BeamParticleData.h"
#include "Herwig/Shower/ShowerAlpha.h"
#include "Herwig/PDT/StandardMatchers.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "HwRemDecayer.fh"
namespace Herwig {
using namespace ThePEG;
/**
* The HwRemDecayer class is responsible for the decay of the remnants. Additional
* secondary scatters have to be evolved backwards to a gluon, the
* first/hard interaction has to be evolved back to a valence quark.
* This is all generated inside this class,
* which main methods are then called by the ShowerHandler.
*
* A simple forced splitting algorithm is used.
* This takes the Remnant object produced from the PDF and backward
* evolution (hadron - parton) and produce partons with the remaining
* flavours and with the correct colour connections.
*
* The algorithim operates by starting with the parton which enters the hard process.
* If this is from the sea there is a forced branching to produce the antiparticle
* from a gluon branching. If the parton entering the hard process was a gluon, or
* a gluon was produced from the first step of the algorithm, there is then a further
* branching back to a valence parton. After these partons have been produced a quark or
* diquark is produced to give the remaining valence content of the incoming hadron.
*
* The forced branching are generated using a scale between QSpac and EmissionRange times
* the minimum scale. The energy fractions are then distributed using
* \f[\frac{\alpha_S}{2\pi}\frac{P(z)}{z}f(x/z,\tilde{q})\f]
* with the massless splitting functions.
*
* \author Manuel B\"ahr
*
* @see \ref HwRemDecayerInterfaces "The interfaces"
* defined for HwRemDecayer.
*/
class HwRemDecayer: public RemnantDecayer {
public:
/** Typedef to store information about colour partners */
typedef vector<pair<tPPtr, tPPtr> > PartnerMap;
public:
/**
* The default constructor.
*/
HwRemDecayer() : allowTop_(false), multiPeriph_(true), quarkPair_(false),
ptmin_(-1.*GeV), beta_(ZERO),
maxtrySoft_(10),
colourDisrupt_(1.0),
ladderbFactor_(0.0),
ladderPower_(-0.08),
ladderNorm_(1.0),
ladderMult_(1.0),
gaussWidth_(0.1),
valOfN_(0),
initTotRap_(0),
_kinCutoff(0.75*GeV),
_forcedSplitScale(2.5*GeV),
_range(1.1), _zbin(0.05),_ybin(0.),
_nbinmax(100), DISRemnantOpt_(0),
PtDistribution_(0),
pomeronStructure_(0), mg_(ZERO) {}
/** @name Virtual functions required by the Decayer class. */
//@{
/**
* Check if this decayer can perfom the decay specified by the
* given decay mode.
* @return true if this decayer can handle the given mode, otherwise false.
*/
virtual bool accept(const DecayMode &) const {
return true;
}
/**
* Return true if this decayer can handle the extraction of the \a
* extracted parton from the given \a particle.
*/
virtual bool canHandle(tcPDPtr particle, tcPDPtr parton) const;
/**
* Return true if this decayed can extract more than one parton from
* a particle.
*/
virtual bool multiCapable() const {
return true;
}
/**
* Perform a decay for a given DecayMode and a given Particle instance.
* @param dm the DecayMode describing the decay.
* @param p the Particle instance to be decayed.
* @param step the step we are working on.
* @return a ParticleVector containing the decay products.
*/
virtual ParticleVector decay(const DecayMode & dm, const Particle & p, Step & step) const;
//@}
public:
/**
* struct that is used to catch exceptions which are thrown
* due to energy conservation issues of additional soft scatters
*/
struct ExtraSoftScatterVeto {};
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
/**
* Do several checks and initialization, for remnantdecay inside ShowerHandler.
*/
void initialize(pair<tRemPPtr, tRemPPtr> rems, tPPair beam, Step & step,
Energy forcedSplitScale);
/**
* Initialize the soft scattering machinery.
* @param ptmin = the pt cutoff used in the UE model
* @param beta = slope of the soft pt-spectrum
*/
void initSoftInteractions(Energy ptmin, InvEnergy2 beta);
/**
* Perform the acual forced splitting.
* @param partons is a pair of ThePEG::Particle pointers which store the final
* partons on which the shower ends.
* @param pdfs are pointers to the pdf objects for both beams
* @param first is a flage wether or not this is the first or a secondary interation
*/
void doSplit(pair<tPPtr, tPPtr> partons, pair<tcPDFPtr, tcPDFPtr> pdfs, bool first);
/**
* Perform the final creation of the diquarks. Set the remnant masses and do
* all colour connections.
* @param colourDisrupt = variable to control how many "hard" scatters
* are colour isolated
* @param softInt = parameter for the number of soft scatters
*/
void finalize(double colourDisrupt=0.0, unsigned int softInt=0);
/**
* Find the children
*/
void findChildren(tPPtr,vector<PPtr> &) const;
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 an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit() {
Interfaced::doinit();
_ybin=0.25/_zbin;
mg_ = getParticleData(ParticleID::g)->constituentMass();
}
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- HwRemDecayer & operator=(const HwRemDecayer &);
+ HwRemDecayer & operator=(const HwRemDecayer &) = delete;
public:
/**
* Simple struct to store info about baryon quark and di-quark
* constituents.
*/
struct HadronContent {
/**
* manually extract the valence flavour \a id.
*/
inline void extract(int id) {
for(unsigned int i=0; i<flav.size(); i++) {
if(id == sign*flav[i]){
if(hadron->id() == ParticleID::gamma ||
(hadron->id() == ParticleID::pomeron && pomeronStructure==1) ||
hadron->id() == ParticleID::reggeon) {
flav[0] = id;
flav[1] = -id;
extracted = 0;
flav.resize(2);
}
else if (hadron->id() == ParticleID::pomeron && pomeronStructure==0) {
extracted = 0;
}
else {
extracted = i;
}
break;
}
}
}
/**
* Return a proper particle ID assuming that \a id has been removed
* from the hadron.
*/
long RemID() const;
/**
* Method to determine whether \a parton is a quark from the sea.
* @return TRUE if \a parton is neither a valence quark nor a gluon.
*/
bool isSeaQuark(tcPPtr parton) const {
return ((parton->id() != ParticleID::g) && ( !isValenceQuark(parton) ) );
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuark(tcPPtr parton) const {
return isValenceQuark(parton->id());
}
/**
* Method to determine whether \a parton is a quark from the sea.
* @return TRUE if \a parton is neither a valence quark nor a gluon.
*/
bool isSeaQuarkData(tcPDPtr partonData) const {
return ((partonData->id() != ParticleID::g) && ( !isValenceQuarkData(partonData) ) );
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuarkData(tcPDPtr partonData) const {
int id(sign*partonData->id());
return find(flav.begin(),flav.end(),id) != flav.end();
}
/**
* Method to determine whether \a parton is a valence quark.
*/
bool isValenceQuark(int id) const {
return find(flav.begin(),flav.end(),sign*id) != flav.end();
}
/** The valence flavours of the corresponding baryon. */
vector<int> flav;
/** The array index of the extracted particle. */
int extracted;
/** -1 if the particle is an anti-particle. +1 otherwise. */
int sign;
/** The ParticleData objects of the hadron */
tcPDPtr hadron;
/** Pomeron treatment */
unsigned int pomeronStructure;
};
/**
* Return the hadron content objects for the incoming particles.
*/
const pair<HadronContent, HadronContent>& content() const {
return theContent;
}
/**
* Return a HadronContent struct from a PPtr to a hadron.
*/
HadronContent getHadronContent(tcPPtr hadron) const;
/**
* Set the hadron contents.
*/
void setHadronContent(tPPair beam) {
theContent.first = getHadronContent(beam.first);
theContent.second = getHadronContent(beam.second);
}
private:
/**
* Do the forced Splitting of the Remnant with respect to the
* extracted parton \a parton.
* @param parton = PPtr to the parton going into the subprocess.
* @param content = HadronContent struct to keep track of flavours.
* @param rem = Pointer to the ThePEG::RemnantParticle.
* @param used = Momentum vector to keep track of remaining momenta.
* @param partners = Vector of pairs filled with tPPtr to the particles
* which should be colour connected.
* @param pdf pointer to the PDF Object which is used for this particle
* @param first = Flag for the first interaction.
*/
void split(tPPtr parton, HadronContent & content, tRemPPtr rem,
Lorentz5Momentum & used, PartnerMap & partners, tcPDFPtr pdf, bool first);
/**
* Merge the colour lines of two particles
* @param p1 = Pointer to particle 1
* @param p2 = Pointer to particle 2
* @param anti = flag to indicate, if (anti)colour was extracted as first parton.
*/
void mergeColour(tPPtr p1, tPPtr p2, bool anti) const;
/**
* Set the colour connections.
* @param partners = Object that holds the information which particles to connect.
* @param anti = flag to indicate, if (anti)colour was extracted as first parton.
* @param disrupt parameter for disruption of the colour structure
*/
void fixColours(PartnerMap partners, bool anti, double disrupt) const;
/**
* Set the momenta of the Remnants properly and boost the decay particles.
*/
void setRemMasses() const;
/**
* This creates a parton from the remaining flavours of the hadron. The
* last parton used was a valance parton, so only 2 (or 1, if meson) flavours
* remain to be used.
*/
PPtr finalSplit(const tRemPPtr rem, long remID,
Lorentz5Momentum usedMomentum) const {
// Create the remnant and set its momentum, also reset all of the decay
// products from the hadron
PPtr remnant = new_ptr(Particle(getParticleData(remID)));
Lorentz5Momentum prem(rem->momentum()-usedMomentum);
prem.setMass(getParticleData(remID)->constituentMass());
prem.rescaleEnergy();
remnant->set5Momentum(prem);
// Add the remnant to the step, but don't do colour connections
thestep->addDecayProduct(rem,remnant,false);
return remnant;
}
/**
* This takes the particle and find a splitting for np -> p + child and
* creates the correct kinematics and connects for such a split. This
* Splitting has an upper bound on qtilde given by the energy argument
* @param rem The Remnant
* @param child The PDG code for the outgoing particle
* @param oldQ The maximum scale for the evolution
* @param oldx The fraction of the hadron's momentum carried by the last parton
* @param pf The momentum of the last parton at input and after branching at output
* @param p The total emitted momentum
* @param content The content of the hadron
*/
PPtr forceSplit(const tRemPPtr rem, long child, Energy &oldQ, double &oldx,
Lorentz5Momentum &pf, Lorentz5Momentum &p,
HadronContent & content) const;
/**
* Check if a particle is a parton from a hadron or not
* @param parton The parton to be tested
*/
bool isPartonic(tPPtr parton) const;
/** @name Soft interaction methods. */
//@{
/**
* Produce pt values according to dN/dp_T = N p_T exp(-beta_*p_T^2)
*/
Energy softPt() const;
/**
* Get the 2 pairs of 5Momenta for the scattering. Needs calling of
* initSoftInteractions.
*/
void softKinematics(Lorentz5Momentum &r1, Lorentz5Momentum &r2,
Lorentz5Momentum &g1, Lorentz5Momentum &g2) const;
/**
* Create N soft gluon interactions
*/
void doSoftInteractions(unsigned int N){
if(!multiPeriph_){
doSoftInteractions_old(N);} //outdated model for soft interactions
else{
doSoftInteractions_multiPeriph(N); // Multiperipheral model
}
}
/**
* Create N soft gluon interactions (old version)
*/
void doSoftInteractions_old(unsigned int N);
/**
* Create N soft gluon interactions with multiperhpheral kinematics
*/
void doSoftInteractions_multiPeriph(unsigned int N);
/**
* Phase space generation for the ladder partons
*/
bool doPhaseSpaceGenerationGluons(vector<Lorentz5Momentum> &softGluons, Energy energy, unsigned int &its)
const;
/**
* This returns the rotation matrix needed to rotate p into the z axis
*/
LorentzRotation rotate(const LorentzMomentum &p) const;
/**
* Methods to generate random distributions also all stolen form UA5Handler
**/
template <typename T>
inline T gaussDistribution(T mean, T stdev) const{
double x = rnd();
x = sqrt(-2.*log(x));
double y;
randAzm(x,x,y);
return mean + stdev*x;
}
/**
* This returns a random number with a flat distribution
* [-A,A] plus gaussian tail with stdev B
* TODO: Should move this to Utilities
* @param A The width of the flat part
* @param B The standard deviation of the gaussian tail
* @return the randomly generated value
*/
inline double randUng(double A, double B) const{
double prun;
if(A == 0.) prun = 0.;
else prun = 1./(1.+B*1.2533/A);
if(rnd() < prun) return 2.*(rnd()-0.5)*A;
else {
double temp = gaussDistribution(0.,B);
if(temp < 0) return temp - abs(A);
else return temp + abs(A);
}
}
template <typename T>
inline void randAzm(T pt, T &px, T &py) const{
double c,s,cs;
while(true) {
c = 2.*rnd()-1.;
s = 2.*rnd()-1.;
cs = c*c+s*s;
if(cs <= 1.&&cs!=0.) break;
}
T qt = pt/cs;
px = (c*c-s*s)*qt;
py = 2.*c*s*qt;
}
inline Energy randExt(Energy AM0,InvEnergy B) const{
double r = rnd();
// Starting value
Energy am = AM0-log(r)/B;
for(int i = 1; i<20; ++i) {
double a = exp(-B*(am-AM0))/(1.+B*AM0);
double f = (1.+B*am)*a-r;
InvEnergy df = -B*B*am*a;
Energy dam = -f/df;
am += dam;
if(am<AM0) am = AM0 + .001*MeV;
if(abs(dam) < .001*MeV) break;
}
return am;
}
/**
* Method to add a particle to the step
* @param parent = pointer to the parent particle
* @param id = Particle ID of the newly created particle
* @param p = Lorentz5Momentum of the new particle
*/
tPPtr addParticle(tcPPtr parent, long id, Lorentz5Momentum p) const;
//@}
/**
* A flag which indicates, whether the extracted valence quark was a
* anti particle.
*/
pair<bool, bool> theanti;
/**
* variable to sum up the x values of the extracted particles
*/
pair<double, double> theX;
/**Pair of HadronContent structs to know about the quark content of the beams*/
pair<HadronContent, HadronContent> theContent;
/**Pair of Lorentz5Momentum to keep track of the forced splitting product momenta*/
pair<Lorentz5Momentum, Lorentz5Momentum> theUsed;
/**
* Pair of PartnerMap's to store the particles, which will be colour
* connected in the end.
*/
pair<PartnerMap, PartnerMap> theMaps;
/**
* Variable to hold a pointer to the current step. The variable is used to
* determine, wether decay(const DecayMode & dm, const Particle & p, Step & step)
* has been called in this event or not.
*/
StepPtr thestep;
/**
* Pair of Remnant pointers. This is needed to boost
* in the Remnant-Remnant CMF after all have been decayed.
*/
pair<RemPPtr, RemPPtr> theRems;
/**
* The beam particle data for the current incoming hadron
*/
mutable tcPPtr theBeam;
/**
* the beam data
*/
mutable Ptr<BeamParticleData>::const_pointer theBeamData;
/**
* The PDF for the current initial-state shower
*/
mutable tcPDFPtr _pdf;
private:
/**
* Switch to control handling of top quarks in proton
*/
bool allowTop_;
/**
* Switch to control using multiperipheral kinemaics
*/
bool multiPeriph_;
/**
* True if kinematics is to be calculated for quarks
*/
bool quarkPair_;
/** @name Soft interaction variables. */
//@{
/**
* Pair of soft Remnant pointers, i.e. Diquarks.
*/
tPPair softRems_;
/**
* ptcut of the UE model
*/
Energy ptmin_;
/**
* slope of the soft pt-spectrum: dN/dp_T = N p_T exp(-beta*p_T^2)
*/
InvEnergy2 beta_;
/**
* Maximum number of attempts for the regeneration of an additional
* soft scattering, before the number of scatters is reduced.
*/
unsigned int maxtrySoft_;
/**
* Variable to store the relative number of colour disrupted
* connections to additional soft subprocesses.
*/
double colourDisrupt_;
/**
* Variable to store the additive factor of the
multiperipheral ladder multiplicity.
*/
double ladderbFactor_;
/**
* Variable of the parameterization of the ladder multiplicity.
*/
double ladderPower_;
/**
* Variable of the parameterization of the ladder multiplicity.
*/
double ladderNorm_;
double ladderMult_;
/**
* Variable to store the gaussian width of the
* fluctuation of the longitudinal momentum
* fraction.
*/
double gaussWidth_;
/**
* Variable to store the current total multiplicity
of a ladder.
*/
double valOfN_;
/**
* Variable to store the initial total rapidity between
of the remnants.
*/
double initTotRap_;
//@}
/** @name Forced splitting variables. */
//@{
/**
* The kinematic cut-off
*/
Energy _kinCutoff;
/**
* The PDF freezing scale as set in ShowerHandler
*/
Energy _forcedSplitScale;
/**
* Range for emission
*/
double _range;
/**
* Size of the bins in z for the interpolation
*/
double _zbin;
/**
* Size of the bins in y for the interpolation
*/
double _ybin;
/**
* Maximum number of bins for the z interpolation
*/
int _nbinmax;
/**
* Pointer to the object calculating the QCD coupling
*/
ShowerAlphaPtr _alphaS;
/**
* Pointer to the object calculating the QED coupling
*/
ShowerAlphaPtr _alphaEM;
/**
* Option for the DIS remnant
*/
unsigned int DISRemnantOpt_;
/**
* Option for the pT generation
*/
unsigned int PtDistribution_;
/**
* Option for the treatment of the pomeron structure
*/
unsigned int pomeronStructure_;
//@}
/**
* The gluon constituent mass.
*/
Energy mg_;
};
}
#endif /* HERWIG_HwRemDecayer_H */
diff --git a/PDF/MultiPartonExtractor.h b/PDF/MultiPartonExtractor.h
--- a/PDF/MultiPartonExtractor.h
+++ b/PDF/MultiPartonExtractor.h
@@ -1,114 +1,114 @@
// -*- C++ -*-
#ifndef Herwig_MultiPartonExtractor_H
#define Herwig_MultiPartonExtractor_H
//
// This is the declaration of the MultiPartonExtractor class.
//
#include "ThePEG/PDF/PartonExtractor.h"
#include <deque>
namespace Herwig {
using namespace ThePEG;
/**
* The MultiPartonExtractor class inherits from the PartonExtractor of ThePEG
* but allows more control over the PDFs used in the case that there are multiple
* stages of parton extraction
*
* @see \ref MultiPartonExtractorInterfaces "The interfaces"
* defined for MultiPartonExtractor.
*/
class MultiPartonExtractor: public PartonExtractor {
public:
/**
* The default constructor.
*/
MultiPartonExtractor() {};
/**
* Return a vector of possible pairs of parton bins which can be
* produced within a given maximum total particle-particle
* invariant mass squared, \a maxEnergy sBin.
*/
virtual PartonPairVec getPartons(Energy maxEnergy, const cPDPair &,
const Cuts &) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/**
* Add parton bins to pbins for the given incoming particle and the
* specified cuts.
*/
virtual void addPartons(tPBPtr incoming ,const PDFCuts & cuts,
std::deque<tcPDFPtr> pdf ,PartonVector & pbins) const;
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- MultiPartonExtractor & operator=(const MultiPartonExtractor &);
+ MultiPartonExtractor & operator=(const MultiPartonExtractor &) = delete;
/**
* PDFBase object to override first PDF
*/
vector<PDFPtr> firstPDF_;
/**
* PDFBase object to override second PDF
*/
vector<PDFPtr> secondPDF_;
};
}
#endif /* Herwig_MultiPartonExtractor_H */
diff --git a/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h b/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h
--- a/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h
+++ b/Shower/Dipole/Kernels/ColourMatrixElementCorrection.h
@@ -1,171 +1,171 @@
// -*- C++ -*-
//
// ColourMatrixElementCorrection.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef Herwig_ColourMatrixElementCorrection_H
#define Herwig_ColourMatrixElementCorrection_H
//
// This is the declaration of the ColourMatrixElementCorrection class.
//
#include "Herwig/Shower/Dipole/Base/DipoleSplittingReweight.h"
#include <tuple>
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Johan Thoren, Simon Platzer
*
* \brief ColourMatrixElementCorrection is implementing colour matrix element
* corrections through the weighted Sudakov algorithm
*
* @see \ref ColourMatrixElementCorrectionInterfaces "The interfaces"
* defined for ColourMatrixElementCorrection.
*/
class ColourMatrixElementCorrection: public DipoleSplittingReweight {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
ColourMatrixElementCorrection();
/**
* The destructor.
*/
virtual ~ColourMatrixElementCorrection();
//@}
public:
/**
* Calculate and cache the colour matrix element correction factor
* for the given splitting type.
*/
double cmec(const DipoleSplittingInfo&) const;
/**
* Return the reweighting factor for the given splitting type.
*/
virtual double evaluate(const DipoleSplittingInfo& s) const {
return cmec(s);
}
/**
* Return the absolute value of the colour matrix element correction
* as an enhancement hint for the sampling of the un-reweighted
* splitting kernel.
*/
virtual double hint(const DipoleSplittingInfo& s) const {
if ( hintOnly(s) )
return cmec(s);
return abs(cmec(s))*lambda;
}
/**
* Return true, if the reweight can be entirely absorbed into the hint. A
* possible detuning will be switched off.
*/
virtual bool hintOnly(const DipoleSplittingInfo& s) const {
return cmec(s) > 0.;
}
/**
* Set the factor in front of enhance used by the veto algorithm.
*/
virtual void reweightFactor(const double c) {
assert(c > 0.0);
lambda = c;
}
/**
* Set the factor in front of enhance used by the veto algorithm.
*/
virtual void negativeScaling(const double c) {
assert(c >= 0.0);
negCMECScaling = c;
}
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* Factor to shuffle the magnitude of the CMEC between the splitting kernel
* and the weight in the reweighted veto algorithm.
*/
double lambda;
/**
* Scaling factor multiplying all of the negative colour matrix element
* corrections. The physically sensible value is 1.0, but this factor can
* be used to examine the effects of the negative contributions.
*/
double negCMECScaling;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- ColourMatrixElementCorrection & operator=(const ColourMatrixElementCorrection &);
+ ColourMatrixElementCorrection & operator=(const ColourMatrixElementCorrection &) = delete;
};
}
#endif /* Herwig_ColourMatrixElementCorrection_H */
diff --git a/Shower/Dipole/Kinematics/IFMassiveKinematics.h b/Shower/Dipole/Kinematics/IFMassiveKinematics.h
--- a/Shower/Dipole/Kinematics/IFMassiveKinematics.h
+++ b/Shower/Dipole/Kinematics/IFMassiveKinematics.h
@@ -1,279 +1,279 @@
// -*- C++ -*-
//
// IFLightKinematics.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 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_IFLightKinematics_H
#define HERWIG_IFLightKinematics_H
//
// This is the declaration of the IFLightKinematics class.
//
#include "DipoleSplittingKinematics.h"
namespace Herwig {
using namespace ThePEG;
/**
* \ingroup DipoleShower
* \author Simon Platzer, Martin Stoll
*
* \brief IFMassiveKinematics implements massless splittings
* off an initial-final dipole.
*
* @see \ref IFMassiveKinematicsInterfaces "The interfaces"
* defined for IFMassiveKinematics.
*/
class IFMassiveKinematics: public DipoleSplittingKinematics {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
IFMassiveKinematics();
/**
* The destructor.
*/
virtual ~IFMassiveKinematics();
//@}
public:
/**
* Return the boundaries in between the evolution
* variable random number is to be sampled; the lower
* cuoff is assumed to correspond to the infrared cutoff.
*/
virtual pair<double,double> kappaSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries in between the momentum
* fraction random number is to be sampled.
*/
virtual pair<double,double> xiSupport(const DipoleSplittingInfo& dIndex) const;
/**
* Return the boundaries on the momentum fraction
*/
virtual pair<double,double> zBoundaries(Energy,
const DipoleSplittingInfo&,
const DipoleSplittingKernel&) const {
return {0.0,1.0};
}
/**
* Return the dipole scale associated to the
* given pair of emitter and spectator. This
* should be the invariant mass or absolute value
* final/final or initial/initial and the absolute
* value of the momentum transfer for intial/final or
* final/initial dipoles.
*/
virtual Energy dipoleScale(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double emX, double specX,
const DipoleSplittingInfo& dInfo,
const DipoleSplittingKernel& split) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy dScale,
double, double,
const DipoleIndex& dIndex,
const DipoleSplittingKernel& split,
tPPtr emitter, tPPtr) const;
/**
* Return the maximum pt for the given dipole scale.
*/
virtual Energy ptMax(Energy,
double, double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
// Only the DipoleSplittingInfo version should be used for massive
// dipoles, for now anyway.
assert(false);
return ZERO;
}
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy dScale,
double emX, double specX,
const DipoleSplittingInfo& dInfo,
const DipoleSplittingKernel& split) const;
/**
* Return the maximum virtuality for the given dipole scale.
*/
virtual Energy QMax(Energy,
double, double,
const DipoleIndex&,
const DipoleSplittingKernel&) const {
// Only the DipoleSplittingInfo version should be used for massive
// dipoles, for now anyway.
assert(false);
return ZERO;
}
/**
* Return the pt given a virtuality.
*/
virtual Energy PtFromQ(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the virtuality given a pt.
*/
virtual Energy QFromPt(Energy scale, const DipoleSplittingInfo&) const;
/**
* Return the random number associated to
* the given pt.
*/
virtual double ptToRandom(Energy pt, Energy dScale,
double emX, double specX,
const DipoleIndex& dIndex,
const DipoleSplittingKernel&) const;
/**
* Generate splitting variables given three random numbers
* and the momentum fractions of the emitter and spectator.
* Return true on success.
*/
virtual bool generateSplitting(double kappa, double xi, double phi,
DipoleSplittingInfo& dIndex,
const DipoleSplittingKernel&);
/**
* Generate the full kinematics given emitter and
* spectator momentum and a previously completeted
* DipoleSplittingInfo object.
*/
virtual void generateKinematics(const Lorentz5Momentum& pEmitter,
const Lorentz5Momentum& pSpectator,
const DipoleSplittingInfo& dInfo);
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** 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;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<IFMassiveKinematics> initIFMassiveKinematics;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- IFMassiveKinematics & operator=(const IFMassiveKinematics &);
+ IFMassiveKinematics & operator=(const IFMassiveKinematics &) = delete;
private:
/**
* Wether or not to choose the `collinear' scheme
*/
bool theCollinearScheme;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of IFMassiveKinematics. */
template <>
struct BaseClassTrait<Herwig::IFMassiveKinematics,1> {
/** Typedef of the first base class of IFMassiveKinematics. */
typedef Herwig::DipoleSplittingKinematics NthBase;
};
/** This template specialization informs ThePEG about the name of
* the IFMassiveKinematics class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::IFMassiveKinematics>
: public ClassTraitsBase<Herwig::IFMassiveKinematics> {
/** Return a platform-independent class name */
static string className() { return "Herwig::IFMassiveKinematics"; }
/**
* The name of a file containing the dynamic library where the class
* IFMassiveKinematics is implemented. It may also include several, space-separated,
* libraries if the class IFMassiveKinematics depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_IFMassiveKinematics_H */
diff --git a/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h b/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h
--- a/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h
+++ b/Shower/Dipole/SpinCorrelations/DipoleShowerParticle.h
@@ -1,182 +1,182 @@
// -*- C++ -*-
//
// DipoleShowerParticle.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleShowerParticle_H
#define HERWIG_DipoleShowerParticle_H
//
// This is the declaration of the DipoleShowerParticle class.
//
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Helicity/WaveFunction/WaveFunctionBase.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "DipoleShowerVertex.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup DipoleShower
*
* \author Stephen Webster
*
*/
class DipoleShowerParticle : public Base {
public:
/**
* Default constructor
*/
DipoleShowerParticle() {}
/**
* Default destructor
**/
~DipoleShowerParticle() {}
public:
/**
* Reset the member variables of the object.
*/
void clear();
/**
* Set up the decay vertex for the emitter.
*/
void prepare( PPtr& part,
const Helicity::Direction emmDir,
const Helicity::Direction specDir,
const Lorentz5Momentum& pVector,
const Lorentz5Momentum& nVector );
/**
* Return the associated decay vertex,
* a DipoleShowerVertex.
**/
DSVertexPtr decayVertex() { return theDecayVertex; }
/**
* Create fermion decay basis states.
* It returns the decay basis states in the
* decay frame as required for the mapping.
*/
vector<LorentzSpinor<SqrtEnergy> > createFermionDecayStates();
/**
* Create vector decay basis states.
*/
void createVectorDecayStates();
/**
* Create fermion production basis states
* for the given particle produced in the splitting.
*/
void createNewFermionSpinInfo( PPtr& outgoing, Helicity::Direction dir);
/**
* Create vector production basis states
* for the given particle produced in the splitting.
*/
void createNewVectorSpinInfo( PPtr& outgoing, Helicity::Direction dir);
/**
* Create the mappings between the production
* and decay states for the fermion and
* store them in the associated decay vertex.
* (No longer applicable) reason for passing the
* decay states as an argument:
* Previously used a check on zero values for computing
* the mapping, rather than a </> 1e-5, this would only
* work when using the original decay state as calculated
* in the decay frame (i.e. without transforming to the
* lab frame and back). Now it simply avoids doing an
* unnecessary rotation of the decay basis
*/
void setFermionMapping( const vector<LorentzSpinor<SqrtEnergy>>& decayBasis );
/**
* Create the mappings between the production
* and decay states the boson and
* store them in the associated decay vertex.
*/
void setVectorMapping();
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
private:
/**
* The pptr to this particle.
*/
PPtr theParticle;
/**
* The dipole shower vertex associated
* with this particle.
*/
DSVertexPtr theDecayVertex;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- DipoleShowerParticle & operator=(const DipoleShowerParticle &);
+ DipoleShowerParticle & operator=(const DipoleShowerParticle &) = delete;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleShowerParticle. */
template <>
struct BaseClassTrait<Herwig::DipoleShowerParticle,1> {
/** Typedef of the first base class of DipoleShowerParticle. */
typedef Base NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleShowerParticle class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleShowerParticle>
: public ClassTraitsBase<Herwig::DipoleShowerParticle> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleShowerParticle"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleShowerParticle is implemented. It may also include several, space-separated,
* libraries if the class DipoleShowerParticle depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleShowerParticle_H */
diff --git a/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h b/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h
--- a/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h
+++ b/Shower/Dipole/SpinCorrelations/DipoleShowerVertex.h
@@ -1,285 +1,285 @@
// -*- C++ -*-
//
// DipoleShowerVertex.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_DipoleShowerVertex_H
#define HERWIG_DipoleShowerVertex_H
//
// This is the declaration of the DipoleShowerVertex class.
//
#include "ThePEG/EventRecord/HelicityVertex.h"
#include "Herwig/Decay/DecayMatrixElement.h"
#include "DipoleShowerVertex.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup DipoleShower
*
* This class represents the vertex for a given splitting
* in the dipole shower.
*
* \author Stephen Webster
*
*/
class DipoleShowerVertex: public HelicityVertex {
public:
/**
* Default constructor
*/
DipoleShowerVertex();
/**
* Default destructor
**/
~DipoleShowerVertex() {}
public:
/**
* Return the matrix element for this vertex.
*/
inline const DecayMEPtr ME() const {
return theMatrixElement;
}
/**
* Set the matrix element
*/
inline void ME(DecayMEPtr in) {
theMatrixElement = in;
}
public:
/**
* Method to calculate the \f$\rho\f$ matrix for one of the decay products
* in the frame of this splitting vertex.
*
* @param iprod The splitting product to compute the \f$\rho\f$ matrix for.
*/
RhoDMatrix getRhoMatrix(int iprod, bool ) const;
/**
* Method to calculate the \f$D\f$ matrix for the decaying
* particle / the incoming to the vertex, in the frame of
* the vertex. The argument is a dummy argument.
*/
RhoDMatrix getDMatrix(int) const;
/**
* Get the lorentz rotation from the working frame
* to the frame of the splitting.
*/
LorentzRotation boostToSplitting();
/**
* Set the p vector for this splitting
*/
void pVector( const Lorentz5Momentum& emitterMom ) { thePVector = emitterMom; }
/**
* Set the n vector for this splitting
*/
void nVector( const Lorentz5Momentum& nMom ) { theNVector = nMom; }
/**
* Set the emitter,Spectator Config (II,IF,FF,FI - F=true, I=false)
*/
void dipoleConfig(const pair<bool,bool>& newConfig) { theDipoleConfig = newConfig; }
/**
* Return the p vector for this splitting
*/
Lorentz5Momentum pVector() const { return thePVector; }
/**
* Return the n/spectator vector for this splitting
*/
Lorentz5Momentum nVector() const { return theNVector; }
/**
* Return the emitter,Spectator Config (II,IF,FF,FI - F=true, I=false)
*/
const pair<bool,bool>& dipoleConfig() const { return theDipoleConfig; }
/**
* Set the decay state to production state
* mapping for this vertex.
*/
void mappingD2P( RhoDMatrix& mapping ) { theMappingDecay2Prod = mapping; }
/**
* Return the mapping from the decay
* states to the production states.
*/
RhoDMatrix mappingD2P() { return theMappingDecay2Prod; }
/**
* Set the production state to decay state
* mapping for this vertex.
*/
void mappingP2D( RhoDMatrix& mapping ) { theMappingProd2Decay = mapping; }
/**
* Return the mapping from the production
* states to the decay states.
*/
RhoDMatrix mappingP2D() { return theMappingProd2Decay; }
/**
* Set the new to old spectator mapping
* for this vertex.
*/
void mappingSpecNewToOld( RhoDMatrix& mapping ) { theMappingSpectatorNewToOld = mapping; }
/**
* Return the new to old spectator mapping
* for this vertex.
*/
RhoDMatrix mappingSpecNewToOld() { return theMappingSpectatorNewToOld; }
/**
* Set the new to old spectator mapping
* for this vertex.
*/
void mappingSpecOldToNew( RhoDMatrix& mapping ) { theMappingSpectatorOldToNew = mapping; }
/**
* Return the new to old spectator mapping
* for this vertex.
*/
RhoDMatrix mappingSpecOldToNew() { return theMappingSpectatorOldToNew; }
public:
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
private:
/**
* Storage of the decay matrix element.
*/
DecayMEPtr theMatrixElement;
/**
* The p vector of the 'splitting basis'
* associated with this vertex.
**/
Lorentz5Momentum thePVector;
/**
* The n vector of the 'splitting basis'
* associated with this vertex.
**/
Lorentz5Momentum theNVector;
/**
* Initial/final config {emitter, spectator}
*/
pair<bool,bool> theDipoleConfig;
/**
* An indicator flag to record if the
* boost to shower for this vertex has been done
*/
bool theBoostCalculated;
/**
* The lorentz transformation from the
* working frame to this splitting.
*/
LorentzRotation theBoostToSplitting;
/**
* The mapping from the decay basis states
* to the production basis states.
*/
RhoDMatrix theMappingDecay2Prod;
/**
* The mapping from the production basis states
* to the decay basis states.
*/
RhoDMatrix theMappingProd2Decay;
/**
* The mapping from the new spectator basis
* states to the old spectator basis states.
*/
RhoDMatrix theMappingSpectatorNewToOld;
/**
* The mapping from the old spectator basis
* states to the new spectator basis states.
*/
RhoDMatrix theMappingSpectatorOldToNew;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- DipoleShowerVertex & operator=(const DipoleShowerVertex &);
+ DipoleShowerVertex & operator=(const DipoleShowerVertex &) = delete;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleShowerVertex. */
template <>
struct BaseClassTrait<Herwig::DipoleShowerVertex,1> {
/** Typedef of the first base class of DipoleShowerVertex. */
typedef HelicityVertex NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleShowerVertex class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleShowerVertex>
: public ClassTraitsBase<Herwig::DipoleShowerVertex> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleShowerVertex"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleShowerVertex is implemented. It may also include several, space-separated,
* libraries if the class DipoleShowerVertex depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* HERWIG_DipoleShowerVertex_H */
diff --git a/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h b/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h
--- a/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h
+++ b/Shower/Dipole/SpinCorrelations/DipoleVertexRecord.h
@@ -1,207 +1,207 @@
// -*- C++ -*-
#ifndef Herwig_DipoleVertexRecord_H
#define Herwig_DipoleVertexRecord_H
//
// This is the declaration of the DipoleVertexRecord class.
//
#include "ThePEG/Config/ThePEG.h"
#include "Herwig/Shower/Dipole/Base/DipoleSplittingInfo.h"
#include "Herwig/Shower/Dipole/Base/Dipole.h"
#include "ThePEG/EventRecord/RhoDMatrix.h"
#include "DipoleShowerParticle.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the DipoleVertexRecord class.
*/
class DipoleVertexRecord: public Base {
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
DipoleVertexRecord() {}
/**
* The destructor.
*/
virtual ~DipoleVertexRecord() { clear(); }
//@}
public:
/**
* Prepare the emitter and spectator
* for the spin correlations computations.
**/
void prepareSplitting( const DipoleSplittingInfo& dInfo, const Dipole& dip);
/**
* Correctly initialise the decay matrix
* to a delta matrix for an external particle.
*/
void initDecayMatrix(PPtr& particle, Helicity::Direction dir);
/**
* Compute the spin density matrix for the given emitter.
* This tracks the path between the given emitter and
* the previous emitter, calculating a rho/decay matrix
* at each vertex as appropriate.
*/
RhoDMatrix emitterDensityMatrix(PPtr emitter);
/**
* Generate the spin-correlated azimuthal angle for a splitting.
*/
void generatePhi(DipoleSplittingInfo& dInfo, Dipole& dip);
/**
* Identify the type of particle and use the appropriate function
* to set up the spin info.
* Required for e.g. MPI
*/
void createSpinInfo(PPtr& part,
const Helicity::Direction& dir);
/**
* Create and set up fermion spin info.
* Required for e.g. MPI
*/
void createFermionSpinInfo(PPtr& part,
const Helicity::Direction& dir);
/**
* Create and set up vector spin info.
* Required for e.g. MPI
*/
void createVectorSpinInfo(PPtr& part,
const Helicity::Direction& dir);
/**
* Update the vertex record following a splitting.
*/
void update(const DipoleSplittingInfo& dInfo);
/**
* For spectators. Set new particle spin info the that of the
* old particle. Update the spin info to include any momentum changes.
*/
void updateSpinInfo( PPtr& oldPart,
PPtr& newPart );
/**
* Set the stopUpdate flag in the spin info of a particle
* incoming to the current decay.
*/
void prepareParticleDecay( const PPtr& parent );
/**
* Update the spin info of the incoming to the decay
* following showering of the decay.
*/
void updateParticleDecay();
/**
* SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests.
* Access the emitter info record.
*/
//map<PPtr,DipoleSplittingInfo> emitterInfoRecord() const {
//return theEmitterInfoRecord;
//}
/**
* SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests.
* Add a splitting to the emitter info record.
*/
// void addToRecord(const DipoleSplittingInfo& dInfo) {
// assert(dInfo.emitter());
// theEmitterInfoRecord[dInfo.emitter()] = dInfo;
// }
/**
* Clear the vertex record: Give up ownership
* on any object involved in the evolution.
*/
virtual void clear();
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
private:
/**
* The current emitter.
*/
DipoleShowerParticle theCurrentEmitter;
/**
* SW 06/02/2019: Required for NearestNeighbourDipoleAnalysis tests.
* Record of the splittings as
* required for the testing analysis.
*/
//map<PPtr, DipoleSplittingInfo> theEmitterInfoRecord;
/**
* The spin info of a particle incoming to the decay
* under consideration.
*/
tcSpinPtr theDecayParentSpinInfo;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
- DipoleVertexRecord & operator=(const DipoleVertexRecord &);
+ DipoleVertexRecord & operator=(const DipoleVertexRecord &) = delete;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of DipoleVertexRecord. */
template <>
struct BaseClassTrait<Herwig::DipoleVertexRecord,1> {
/** Typedef of the first base class of DipoleVertexRecord. */
typedef Base NthBase;
};
/** This template specialization informs ThePEG about the name of
* the DipoleVertexRecord class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::DipoleVertexRecord>
: public ClassTraitsBase<Herwig::DipoleVertexRecord> {
/** Return a platform-independent class name */
static string className() { return "Herwig::DipoleVertexRecord"; }
/**
* The name of a file containing the dynamic library where the class
* DipoleVertexRecord is implemented. It may also include several, space-separated,
* libraries if the class DipoleVertexRecord depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwDipoleShower.so"; }
};
/** @endcond */
}
#endif /* Herwig_DipoleVertexRecord_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 8:34 PM (1 d, 5 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805084
Default Alt Text
(252 KB)

Event Timeline