Page MenuHomeHEPForge
Files F8723342

No OneTemporary
Actions

View Options

diff --git a/Decay/General/SFFDecayer.cc b/Decay/General/SFFDecayer.cc
--- a/Decay/General/SFFDecayer.cc
+++ b/Decay/General/SFFDecayer.cc
@@ -1,493 +1,500 @@
// -*- C++ -*-
//
// SFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SFFDecayer class.
//
#include "SFFDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.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/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "Herwig/Utilities/Kinematics.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
IBPtr SFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr SFFDecayer::fullclone() const {
return new_ptr(*this);
}
void SFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr> vertex,
map<ShowerInteraction,VertexBasePtr> & inV,
const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
map<ShowerInteraction,VertexBasePtr> ) {
decayInfo(incoming,outgoing);
for(auto vert : vertex) {
vertex_ .push_back(dynamic_ptr_cast<AbstractFFSVertexPtr>(vert));
perturbativeVertex_ .push_back(dynamic_ptr_cast<FFSVertexPtr> (vert));
}
vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
for(auto & inter : itemp) {
incomingVertex_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(inV.at(inter));
outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
}
}
void SFFDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << perturbativeVertex_
<< incomingVertex_ << outgoingVertex1_
<< outgoingVertex2_;
}
void SFFDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> perturbativeVertex_
>> incomingVertex_ >> outgoingVertex1_
>> outgoingVertex2_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<SFFDecayer,GeneralTwoBodyDecayer>
+DescribeClass<SFFDecayer,GeneralTwoBodyDecayer2>
describeHerwigSFFDecayer("Herwig::SFFDecayer", "Herwig.so");
void SFFDecayer::Init() {
static ClassDocumentation<SFFDecayer> documentation
("This class implements to decay of a scalar to 2 fermions");
}
-
-double SFFDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,MEOption meopt) const {
- if(!ME())
- ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
- // work out which is the fermion and antifermion
+void SFFDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
int iferm(1),ianti(0);
int itype[2];
for(unsigned int ix=0;ix<2;++ix) {
if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
else itype[ix] = 2;
}
if(itype[0]==0||itype[1]==1||(itype[0]==2&&itype[1]==2)) swap(iferm,ianti);
+ ScalarWaveFunction::
+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),incoming,true);
+ SpinorBarWaveFunction::
+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
+ SpinorWaveFunction::
+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
+}
+double SFFDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ MEOption meopt) const {
+ if(!ME())
+ ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1Half,PDT::Spin1Half)));
+ // work out which is the fermion and antifermion
+ int iferm(1),ianti(0);
+ int itype[2];
+ for(unsigned int ix=0;ix<2;++ix) {
+ if(outgoing[ix]->CC()) itype[ix] = outgoing[ix]->id()>0 ? 0:1;
+ else itype[ix] = 2;
+ }
+ if(itype[0]==0||itype[1]==1||(itype[0]==2&&itype[1]==2)) swap(iferm,ianti);
if(meopt==Initialize) {
ScalarWaveFunction::
- calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),incoming);
- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- ScalarWaveFunction::
- constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
- SpinorBarWaveFunction::
- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
- SpinorWaveFunction::
- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
- return 0.;
- }
SpinorBarWaveFunction::
- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
SpinorWaveFunction::
- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
- Energy2 scale(sqr(inpart.mass()));
+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
+ Energy2 scale(sqr(part.mass()));
for(unsigned int ifm = 0; ifm < 2; ++ifm){
for(unsigned int ia = 0; ia < 2; ++ia) {
if(iferm > ianti) (*ME())(0, ia, ifm) = 0.;
else (*ME())(0, ifm, ia) = 0.;
for(auto vert : vertex_) {
if(iferm > ianti){
(*ME())(0, ia, ifm) += vert->evaluate(scale,wave_[ia],
wavebar_[ifm],swave_);
}
else {
(*ME())(0, ifm, ia) += vert->evaluate(scale,wave_[ia],
wavebar_[ifm],swave_);
}
}
}
}
double output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
// colour and identical particle factors
- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
- decay[1]->dataPtr());
+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
// return the answer
return output;
}
Energy SFFDecayer::partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
if( inpart.second < outa.second + outb.second ) return ZERO;
if(perturbativeVertex_.size()==1 &&
perturbativeVertex_[0]) {
tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outb.first, outa.first,
in);
double mu1(outa.second/inpart.second),mu2(outb.second/inpart.second);
double c2 = norm(perturbativeVertex_[0]->norm());
Complex al(perturbativeVertex_[0]->left()), ar(perturbativeVertex_[0]->right());
double me2 = -c2*( (norm(al) + norm(ar))*( sqr(mu1) + sqr(mu2) - 1.)
+ 2.*(ar*conj(al) + al*conj(ar)).real()*mu1*mu2 );
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
outb.second);
Energy output = me2*pcm/(8*Constants::pi);
// colour factor
output *= colourFactor(inpart.first,outa.first,outb.first);
return output;
}
else {
- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
}
}
double SFFDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
// work out which is the fermion and antifermion
int ianti(0), iferm(1), iglu(2);
int itype[2];
for(unsigned int ix=0;ix<2;++ix) {
if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
else itype[ix] = 2;
}
if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(meopt==Initialize) {
// create scalar wavefunction for decaying particle
ScalarWaveFunction::
calculateWaveFunctions(rho3_,const_ptr_cast<tPPtr>(&inpart),incoming);
swave3_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
}
// setup spin information when needed
if(meopt==Terminate) {
ScalarWaveFunction::
constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorBarWaveFunction::
constructSpinInfo(wavebar3_ ,decay[iferm],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
VectorWaveFunction::
constructSpinInfo(gluon_ ,decay[iglu ],outgoing,true,false);
return 0.;
}
// calculate colour factors and number of colour flows
unsigned int nflow;
vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);
vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, PDT::Spin1Half,
PDT::Spin1Half, PDT::Spin1)));
// create wavefunctions
SpinorBarWaveFunction::
calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
SpinorWaveFunction::
calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
VectorWaveFunction::
calculateWaveFunctions(gluon_ , decay[iglu ],outgoing,true);
// gauge invariance test
#ifdef GAUGE_CHECK
gluon_.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) gluon_.push_back(VectorWaveFunction());
else {
gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
decay[iglu ]->dataPtr(),10,
outgoing));
}
}
#endif
// identify fermion and/or anti-fermion vertex
AbstractFFVVertexPtr outgoingVertexF;
AbstractFFVVertexPtr outgoingVertexA;
identifyVertices(iferm, ianti, inpart, decay, outgoingVertexF, outgoingVertexA,
inter);
- const GeneralTwoBodyDecayer::CFlow & colourFlow
+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
= colourFlows(inpart, decay);
Energy2 scale(sqr(inpart.mass()));
double gs(0.);
bool couplingSet(false);
#ifdef GAUGE_CHECK
double total=0.;
#endif
for(unsigned int ifm = 0; ifm < 2; ++ifm) {
for(unsigned int ia = 0; ia < 2; ++ia) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from the incoming scalar
if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(incomingVertex_[inter]);
ScalarWaveFunction scalarInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),
gluon_[2*ig],swave3_,inpart.mass());
assert(swave3_.particle()->id()==scalarInter.particle()->id());
if(!couplingSet) {
gs = abs(incomingVertex_[inter]->norm());
couplingSet = true;
}
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ia],
wavebar3_[ifm],scalarInter);
for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
(*ME[colourFlow[0][ix].first])(0, ia, ifm, ig) +=
colourFlow[0][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing fermion
if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED)) {
assert(outgoingVertexF);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iferm]->dataPtr();
if(off->CC()) off = off->CC();
SpinorBarWaveFunction interS =
outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
gluon_[2*ig],decay[iferm]->mass());
assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());
if(!couplingSet) {
gs = abs(outgoingVertexF->norm());
couplingSet = true;
}
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ia], interS,swave3_);
for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
(*ME[colourFlow[1][ix].first])(0, ia, ifm, ig) +=
colourFlow[1][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing antifermion
if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED)) {
assert(outgoingVertexA);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[ianti]->dataPtr();
if(off->CC()) off = off->CC();
SpinorWaveFunction interS =
outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
gluon_[2*ig],decay[ianti]->mass());
assert(wave3_[ia].particle()->id()==interS.particle()->id());
if(!couplingSet) {
gs = abs(outgoingVertexA->norm());
couplingSet = true;
}
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,interS,wavebar3_[ifm],swave3_);
for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
(*ME[colourFlow[2][ix].first])(0, ia, ifm, ig) +=
colourFlow[2][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
}
}
}
// contract matrices
double output=0.;
for(unsigned int ix=0; ix<nflow; ++ix){
for(unsigned int iy=0; iy<nflow; ++iy){
output+=cfactors[ix][iy]*(ME[ix]->contract(*ME[iy],rho3_)).real();
}
}
// divide by alpha(S,EM)
output *= (4.*Constants::pi)/sqr(gs);
#ifdef GAUGE_CHECK
double ratio = output/total;
if(abs(ratio)>1e-20) {
generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
for(unsigned int ix=0;ix<decay.size();++ix)
generator()->log() << *decay[ix] << "\n";
generator()->log() << "Test of gauge invariance " << ratio << "\n";
}
#endif
// return the answer
return output;
}
void SFFDecayer::identifyVertices(const int iferm, const int ianti,
const Particle & inpart, const ParticleVector & decay,
AbstractFFVVertexPtr & outgoingVertexF,
AbstractFFVVertexPtr & outgoingVertexA,
ShowerInteraction inter) {
// QCD
if(inter==ShowerInteraction::QCD) {
// work out which fermion each outgoing vertex corresponds to
// two outgoing vertices
if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
((decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) ||
(decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour8))) {
if(outgoingVertex1_[inter]==outgoingVertex2_[inter]) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) {
if(outgoingVertex1_[inter]==outgoingVertex2_[inter]) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
// one outgoing vertex
else if(inpart.dataPtr()->iColour()==PDT::Colour3){
if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour0){
if (outgoingVertex1_[inter]) outgoingVertexF = outgoingVertex1_[inter];
else if(outgoingVertex2_[inter]) outgoingVertexF = outgoingVertex2_[inter];
}
else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour8) {
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr()))) {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
else {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
}
else if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) {
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
}
else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
if(decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour0){
if (outgoingVertex1_[inter]) outgoingVertexA = outgoingVertex1_[inter];
else if(outgoingVertex2_[inter]) outgoingVertexA = outgoingVertex2_[inter];
}
else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
else if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3) {
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
}
else if(inpart.dataPtr()->iColour()==PDT::Colour6 ||
inpart.dataPtr()->iColour()==PDT::Colour6bar) {
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
if (! ((incomingVertex_[inter] && (outgoingVertexF || outgoingVertexA)) ||
( outgoingVertexF && outgoingVertexA))) {
throw Exception()
<< "Invalid vertices for QCD radiation in SFF decay in SFFDecayer::identifyVertices"
<< Exception::runerror;
}
}
// QED
else {
if(decay[iferm]->dataPtr()->charged()) {
if (outgoingVertex1_[inter] &&
outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr())))
outgoingVertexF = outgoingVertex1_[inter];
else
outgoingVertexF = outgoingVertex2_[inter];
}
if(decay[ianti]->dataPtr()->charged()) {
if (outgoingVertex1_[inter] &&
outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr())))
outgoingVertexA = outgoingVertex1_[inter];
else
outgoingVertexA = outgoingVertex2_[inter];
}
}
}
diff --git a/Decay/General/SFFDecayer.h b/Decay/General/SFFDecayer.h
--- a/Decay/General/SFFDecayer.h
+++ b/Decay/General/SFFDecayer.h
@@ -1,234 +1,241 @@
// -*- C++ -*-
//
// SFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SFFDecayer_H
#define HERWIG_SFFDecayer_H
//
// This is the declaration of the SFFDecayer class.
//
-#include "GeneralTwoBodyDecayer.h"
+#include "GeneralTwoBodyDecayer2.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Helicity/Vertex/Scalar/FFSVertex.h"
#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::FFSVertexPtr;
/** \ingroup Decay
* The SFFDecayer class implements the decay of a scalar to 2
* fermions in a general model. It holds an FFSVertex pointer that
* must be typecast from the VertexBase pointer held in
- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
* partialWidth().
*
- * @see GeneralTwoBodyDecayer
+ * @see GeneralTwoBodyDecayer2
*/
-class SFFDecayer: public GeneralTwoBodyDecayer {
+class SFFDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
SFFDecayer() {}
-
- /** @name Virtual functions required by the Decayer class. */
- //@{
- /**
+
+ /**
* 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 ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
- * @param decay The particles produced in the decay.
+ * @param outgoing The particles produced in the decay
+ * @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ 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;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa One of the decay products.
* @param outb The other decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {
POWHEGType output = FSR;
for(auto vertex : vertex_) {
if(vertex->orderInAllCouplings()!=1) {
output = No;
break;
}
}
return output;
}
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter,
MEOption meopt);
/**
* Indentify outgoing vertices for the fermion and antifermion
*/
void identifyVertices(const int iferm, const int ianti,
const Particle & inpart, const ParticleVector & decay,
AbstractFFVVertexPtr & outgoingVertexF,
AbstractFFVVertexPtr & outgoingVertexA,
ShowerInteraction inter);
/**
* Set the information on the decay
*/
virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr>,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > &,
map<ShowerInteraction,VertexBasePtr>);
//@}
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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SFFDecayer & operator=(const SFFDecayer &);
private:
/**
* Abstract pointer to AbstractFFSVertex
*/
vector<AbstractFFSVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<FFSVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractVSSVertex for QCD radiation from incoming scalar
*/
map<ShowerInteraction,AbstractVSSVertexPtr> incomingVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Scalar wavefunction
*/
mutable ScalarWaveFunction swave_;
/**
* Spinor wavefunction
*/
mutable vector<SpinorWaveFunction> wave_;
/**
* Barred spinor wavefunction
*/
mutable vector<SpinorBarWaveFunction> wavebar_;
/**
* Spin density matrix for 3 body decay
*/
mutable RhoDMatrix rho3_;
/**
* Scalar wavefunction for 3 body decay
*/
mutable ScalarWaveFunction swave3_;
/**
* Spinor wavefunction for 3 body decay
*/
mutable vector<SpinorWaveFunction> wave3_;
/**
* Barred spinor wavefunction for 3 body decay
*/
mutable vector<SpinorBarWaveFunction> wavebar3_;
/**
* Vector wavefunction for 3 body decay
*/
mutable vector<VectorWaveFunction> gluon_;
};
}
#endif /* HERWIG_SFFDecayer_H */
diff --git a/Decay/General/SSVDecayer.cc b/Decay/General/SSVDecayer.cc
--- a/Decay/General/SSVDecayer.cc
+++ b/Decay/General/SSVDecayer.cc
@@ -1,368 +1,372 @@
// -*- C++ -*-
//
// SSVDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SSVDecayer class.
//
#include "SSVDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Helicity/WaveFunction/ScalarWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
IBPtr SSVDecayer::clone() const {
return new_ptr(*this);
}
IBPtr SSVDecayer::fullclone() const {
return new_ptr(*this);
}
void SSVDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr> vertex,
map<ShowerInteraction,VertexBasePtr> & inV,
const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
map<ShowerInteraction,VertexBasePtr> fourV) {
decayInfo(incoming,outgoing);
for(auto vert : vertex) {
vertex_ .push_back(dynamic_ptr_cast<AbstractVSSVertexPtr>(vert));
perturbativeVertex_.push_back(dynamic_ptr_cast<VSSVertexPtr> (vert));
}
vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
for(auto & inter : itemp) {
incomingVertex_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(inV.at(inter));
fourPointVertex_[inter] = dynamic_ptr_cast<AbstractVVSSVertexPtr>(fourV.at(inter));
outgoingVertexS_[inter] = AbstractVSSVertexPtr();
outgoingVertexV_[inter] = AbstractVVVVertexPtr();
if(outV[0].at(inter)) {
if (outV[0].at(inter)->getName()==VertexType::VSS)
outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[0].at(inter));
else
outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[0].at(inter));
}
if(outV[1].at(inter)) {
if (outV[1].at(inter)->getName()==VertexType::VSS)
outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
else
outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
}
}
}
void SSVDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << perturbativeVertex_
<< incomingVertex_ << outgoingVertexS_
<< outgoingVertexV_ << fourPointVertex_;
}
void SSVDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> perturbativeVertex_
>> incomingVertex_ >> outgoingVertexS_
>> outgoingVertexV_ >> fourPointVertex_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<SSVDecayer,GeneralTwoBodyDecayer>
+DescribeClass<SSVDecayer,GeneralTwoBodyDecayer2>
describeHerwigSSVDecayer("Herwig::SSVDecayer", "Herwig.so");
void SSVDecayer::Init() {
static ClassDocumentation<SSVDecayer> documentation
("This implements the decay of a scalar to a vector and a scalar");
}
-double SSVDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,
+void SSVDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ unsigned int isc(0),ivec(1);
+ if(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec);
+ ScalarWaveFunction::
+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),incoming,true);
+ ScalarWaveFunction::
+ constructSpinInfo(decay[isc],outgoing,true);
+ VectorWaveFunction::
+ constructSpinInfo(vector_,decay[ivec],outgoing,true,false);
+}
+
+double SSVDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
unsigned int isc(0),ivec(1);
- if(decay[0]->dataPtr()->iSpin() != PDT::Spin0) swap(isc,ivec);
+ if(outgoing[0]->iSpin() != PDT::Spin0) swap(isc,ivec);
if(!ME()) {
if(ivec==1)
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin0,PDT::Spin1)));
else
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin0,PDT::Spin1,PDT::Spin0)));
}
if(meopt==Initialize) {
ScalarWaveFunction::
- calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),incoming);
- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),incoming);
+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),incoming);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- ScalarWaveFunction::
- constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
- ScalarWaveFunction::
- constructSpinInfo(decay[isc],outgoing,true);
- VectorWaveFunction::
- constructSpinInfo(vector_,decay[ivec],outgoing,true,false);
- }
VectorWaveFunction::
- calculateWaveFunctions(vector_,decay[ivec],outgoing,false);
- ScalarWaveFunction sca(decay[isc]->momentum(),decay[isc]->dataPtr(),outgoing);
- Energy2 scale(sqr(inpart.mass()));
+ calculateWaveFunctions(vector_,momenta[ivec],outgoing[ivec],Helicity::outgoing,false);
+ ScalarWaveFunction sca(momenta[isc],outgoing[isc],Helicity::outgoing);
+ Energy2 scale(sqr(part.mass()));
//make sure decay matrix element is in the correct order
double output(0.);
if(ivec == 0) {
for(unsigned int ix = 0; ix < 3; ++ix) {
(*ME())(0, ix, 0) = 0.;
for(auto vert : vertex_)
(*ME())(0, ix, 0) += vert->evaluate(scale,vector_[ix],sca, swave_);
}
}
else {
for(unsigned int ix = 0; ix < 3; ++ix) {
(*ME())(0, 0, ix) = 0.;
for(auto vert : vertex_)
(*ME())(0, 0, ix) += vert->evaluate(scale,vector_[ix],sca,swave_);
}
}
output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
// colour and identical particle factors
- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
- decay[1]->dataPtr());
+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
// return the answer
return output;
}
Energy SSVDecayer:: partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
if( inpart.second < outa.second + outb.second ) return ZERO;
if(perturbativeVertex_.size()==1 &&
perturbativeVertex_[0]) {
double mu1sq(sqr(outa.second/inpart.second)),
mu2sq(sqr(outb.second/inpart.second));
tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
if(outa.first->iSpin() == PDT::Spin0) {
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outb.first, outa.first,in);
}
else {
swap(mu1sq,mu2sq);
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outa.first, outb.first,in);
}
double me2(0.);
if(mu2sq == 0.)
me2 = -2.*mu1sq - 2.;
else
me2 = ( sqr(mu2sq - mu1sq) - 2.*(mu2sq + mu1sq) + 1. )/mu2sq;
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
outb.second);
Energy output = pcm*me2*norm(perturbativeVertex_[0]->norm())/8./Constants::pi;
// colour factor
output *= colourFactor(inpart.first,outa.first,outb.first);
// return the answer
return output;
}
else {
- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
}
}
double SSVDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
int iscal (0), ivect (1), iglu (2);
// get location of outgoing scalar/vector
if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,ivect);
if(meopt==Initialize) {
// create scalar wavefunction for decaying particle
ScalarWaveFunction::calculateWaveFunctions(rho3_,const_ptr_cast<tPPtr>(&inpart),incoming);
swave3_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),incoming);
}
// setup spin information when needed
if(meopt==Terminate) {
ScalarWaveFunction::
constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),incoming,true);
ScalarWaveFunction::
constructSpinInfo(decay[iscal],outgoing,true);
VectorWaveFunction::
constructSpinInfo(vector3_,decay[ivect],outgoing,true,false);
VectorWaveFunction::
constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
return 0.;
}
// calculate colour factors and number of colour flows
unsigned int nflow;
vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);
vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin0, PDT::Spin0,
PDT::Spin1, PDT::Spin1)));
// create wavefunctions
ScalarWaveFunction scal(decay[iscal]->momentum(), decay[iscal]->dataPtr(),outgoing);
VectorWaveFunction::calculateWaveFunctions(vector3_,decay[ivect],outgoing,false);
VectorWaveFunction::calculateWaveFunctions(gluon_, decay[iglu ],outgoing,true );
// gauge invariance test
#ifdef GAUGE_CHECK
gluon_.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) gluon_.push_back(VectorWaveFunction());
else {
gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
decay[iglu ]->dataPtr(),10,
outgoing));
}
}
#endif
if (! ((incomingVertex_[inter] && (outgoingVertexS_[inter] || outgoingVertexV_[inter])) ||
(outgoingVertexS_[inter] && outgoingVertexV_[inter])))
throw Exception()
<< "Invalid vertices for radiation in SSV decay in SSVDecayer::threeBodyME"
<< Exception::runerror;
// sort out colour flows
int S(1), V(2);
if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[ivect]->dataPtr()->iColour()==PDT::Colour8)
swap(S,V);
else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 &&
decay[iscal]->dataPtr()->iColour()==PDT::Colour8)
swap(S,V);
Energy2 scale(sqr(inpart.mass()));
- const GeneralTwoBodyDecayer::CFlow & colourFlow
+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
= colourFlows(inpart, decay);
double gs(0.);
bool couplingSet(false);
#ifdef GAUGE_CHECK
double total=0.;
#endif
for(unsigned int iv = 0; iv < 3; ++iv) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from the incoming scalar
if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(incomingVertex_[inter]);
ScalarWaveFunction scalarInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),
gluon_[2*ig],swave3_,inpart.mass());
assert(swave3_.particle()->id()==scalarInter.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,vector3_[iv],scal,scalarInter);
if(!couplingSet) {
gs = abs(incomingVertex_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
(*ME[colourFlow[0][ix].first])(0, 0, iv, ig) +=
colourFlow[0][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from the outgoing scalar
if((decay[iscal]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[iscal]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexS_[inter]);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iscal]->dataPtr();
if(off->CC()) off = off->CC();
ScalarWaveFunction scalarInter =
outgoingVertexS_[inter]->evaluate(scale,3,off,gluon_[2*ig],scal,decay[iscal]->mass());
assert(scal.particle()->id()==scalarInter.particle()->id());
if(!couplingSet) {
gs = abs(outgoingVertexS_[inter]->norm());
couplingSet = true;
}
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,vector3_[iv],scalarInter,swave3_);
for(unsigned int ix=0;ix<colourFlow[S].size();++ix) {
(*ME[colourFlow[S][ix].first])(0, 0, iv, ig) +=
colourFlow[S][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing vector
if((decay[ivect]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[ivect]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexV_[inter]);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[ivect]->dataPtr();
if(off->CC()) off = off->CC();
VectorWaveFunction vectorInter =
outgoingVertexV_[inter]->evaluate(scale,3,off,gluon_[2*ig],
vector3_[iv],decay[ivect]->mass());
assert(vector3_[iv].particle()->id()==vectorInter.particle()->id());
if(!couplingSet) {
gs = abs(outgoingVertexV_[inter]->norm());
couplingSet = true;
}
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,vectorInter,scal,swave3_);
for(unsigned int ix=0;ix<colourFlow[V].size();++ix) {
(*ME[colourFlow[V][ix].first])(0, 0, iv, ig) +=
colourFlow[V][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from 4 point vertex
if (fourPointVertex_[inter]) {
Complex diag = fourPointVertex_[inter]->evaluate(scale, gluon_[2*ig], vector3_[iv],
scal, swave3_);
for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
(*ME[colourFlow[3][ix].first])(0, 0, iv, ig) +=
colourFlow[3][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
}
}
// contract matrices
double output=0.;
for(unsigned int ix=0; ix<nflow; ++ix){
for(unsigned int iy=0; iy<nflow; ++iy){
output+=cfactors[ix][iy]*(ME[ix]->contract(*ME[iy],rho3_)).real();
}
}
// divide by alpha_(S,EM)
output*=(4.*Constants::pi)/sqr(gs);
#ifdef GAUGE_CHECK
double ratio = output/total;
if(abs(ratio)>1e-20) {
generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
for(unsigned int ix=0;ix<decay.size();++ix)
generator()->log() << *decay[ix] << "\n";
generator()->log() << "Test of gauge invariance " << ratio << "\n";
}
#endif
// return the answer
return output;
}
diff --git a/Decay/General/SSVDecayer.h b/Decay/General/SSVDecayer.h
--- a/Decay/General/SSVDecayer.h
+++ b/Decay/General/SSVDecayer.h
@@ -1,224 +1,233 @@
// -*- C++ -*-
//
// SSVDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_SSVDecayer_H
#define HERWIG_SSVDecayer_H
//
// This is the declaration of the SSVDecayer class.
//
-#include "GeneralTwoBodyDecayer.h"
+#include "GeneralTwoBodyDecayer2.h"
#include "ThePEG/Helicity/Vertex/Scalar/VSSVertex.h"
#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"
#include "ThePEG/Helicity/Vertex/Scalar/VVSSVertex.h"
#include "ThePEG/Repository/EventGenerator.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::VSSVertexPtr;
/** \ingroup Decay
* The SSVDecayer class implements the decay of a scalar to a vector
* and a scalar in a general model. It holds an VSSVertex pointer
* that must be typecast from the VertexBase pointer held in
- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
* partialWidth().
*
- * @see GeneralTwoBodyDecayer
+ * @see GeneralTwoBodyDecayer2
*/
-class SSVDecayer: public GeneralTwoBodyDecayer {
+class SSVDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
SSVDecayer() {}
/** @name Virtual functions required by the Decayer class. */
//@{
/**
- * Return the matrix element squared for a given mode and phase-space channel
- * @param ichan The channel we are calculating the matrix element for.
+ * Return the matrix element squared for a given mode and phase-space channel.
+ * @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
- * @param decay The particles produced in the decay.
+ * @param outgoing The particles produced in the decay
+ * @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ 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;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa One of the decay products.
* @param outb The other decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {
POWHEGType output = FSR;
for(auto vertex : vertex_) {
if(vertex->orderInAllCouplings()!=1) {
output = No;
break;
}
}
return output;
}
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter,
MEOption meopt);
/**
* Set the information on the decay
*/
virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr>,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > &,
map<ShowerInteraction,VertexBasePtr>);
//@}
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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SSVDecayer & operator=(const SSVDecayer &);
private:
/**
* Abstract pointer to AbstractFFVVertex
*/
vector<AbstractVSSVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<VSSVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractVSSVertex for QCD radiation from incoming scalar
*/
map<ShowerInteraction,AbstractVSSVertexPtr> incomingVertex_;
/**
* Abstract pointer to AbstractVSSVertex for QCD radiation from outgoing scalar
*/
map<ShowerInteraction,AbstractVSSVertexPtr> outgoingVertexS_;
/**
* Abstract pointer to AbstractVVVVertex for QCD radiation from outgoing vector
*/
map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertexV_;
/**
* Abstract pointer to AbstractVVSSVertex for QCD radiation from 4 point vertex
*/
map<ShowerInteraction,AbstractVVSSVertexPtr> fourPointVertex_;
/**
* Spinor density matrix
*/
mutable RhoDMatrix rho_;
/**
* Scalar wavefunction
*/
mutable Helicity::ScalarWaveFunction swave_;
/**
* Vector wavefunction
*/
mutable vector<Helicity::VectorWaveFunction> vector_;
/**
* Spin density matrix for 3 body decay
*/
mutable RhoDMatrix rho3_;
/**
* Scalar wavefunction for 3 body decay
*/
mutable Helicity::ScalarWaveFunction swave3_;
/**
* Scalar wavefunction for 3 body decay
*/
mutable Helicity::ScalarWaveFunction scal_;
/**
* Vector wavefunction for 3 body decay
*/
mutable vector<Helicity::VectorWaveFunction> vector3_;
/**
* Vector wavefunction for 3 body decay
*/
mutable vector<Helicity::VectorWaveFunction> gluon_;
};
}
#endif /* HERWIG_SSVDecayer_H */
diff --git a/Decay/General/TFFDecayer.cc b/Decay/General/TFFDecayer.cc
--- a/Decay/General/TFFDecayer.cc
+++ b/Decay/General/TFFDecayer.cc
@@ -1,353 +1,356 @@
// -*- C++ -*-
//
// TFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TFFDecayer class.
//
#include "TFFDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.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/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "Herwig/Utilities/Kinematics.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
IBPtr TFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr TFFDecayer::fullclone() const {
return new_ptr(*this);
}
void TFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr> vertex,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
map<ShowerInteraction,VertexBasePtr> fourV) {
decayInfo(incoming,outgoing);
for(auto vert : vertex) {
vertex_ .push_back(dynamic_ptr_cast<AbstractFFTVertexPtr>(vert));
perturbativeVertex_.push_back(dynamic_ptr_cast<FFTVertexPtr> (vert));
}
vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
for(auto & inter : itemp) {
fourPointVertex_[inter] = dynamic_ptr_cast<AbstractFFVTVertexPtr>(fourV.at(inter));
outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr> (outV[0].at(inter));
outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr> (outV[1].at(inter));
}
}
void TFFDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << perturbativeVertex_
<< outgoingVertex1_ << outgoingVertex2_
<< fourPointVertex_;
}
void TFFDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> perturbativeVertex_
>> outgoingVertex1_ >> outgoingVertex2_
>> fourPointVertex_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<TFFDecayer,GeneralTwoBodyDecayer>
+DescribeClass<TFFDecayer,GeneralTwoBodyDecayer2>
describeHerwigTFFDecayer("Herwig::TFFDecayer", "Herwig.so");
void TFFDecayer::Init() {
static ClassDocumentation<TFFDecayer> documentation
("The TFFDecayer class implements the decay of a tensor particle "
"to 2 fermions ");
}
-double TFFDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,
+void TFFDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ unsigned int iferm(0),ianti(1);
+ if(decay[0]->id()>=0) swap(iferm,ianti);
+ TensorWaveFunction::
+ constructSpinInfo(tensors_,const_ptr_cast<tPPtr>(&part),
+ incoming,true,false);
+ SpinorBarWaveFunction::
+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
+ SpinorWaveFunction::
+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
+}
+
+double TFFDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
unsigned int iferm(0),ianti(1);
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin2,PDT::Spin1Half,PDT::Spin1Half)));
- if(decay[0]->id()>=0) swap(iferm,ianti);
+ if(outgoing[0]->id()>=0) swap(iferm,ianti);
if(meopt==Initialize) {
TensorWaveFunction::
- calculateWaveFunctions(tensors_,rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(tensors_,rho_,const_ptr_cast<tPPtr>(&part),
incoming,false);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- TensorWaveFunction::
- constructSpinInfo(tensors_,const_ptr_cast<tPPtr>(&inpart),
- incoming,true,false);
- SpinorBarWaveFunction::
- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
- SpinorWaveFunction::
- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
- return 0.;
- }
SpinorBarWaveFunction::
- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
SpinorWaveFunction::
- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
- Energy2 scale(sqr(inpart.mass()));
+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
+ Energy2 scale(sqr(part.mass()));
unsigned int thel,fhel,ahel;
for(thel=0;thel<5;++thel) {
for(fhel=0;fhel<2;++fhel) {
for(ahel=0;ahel<2;++ahel) {
if(iferm > ianti) {
(*ME())(thel,fhel,ahel) = 0.;
for(auto vert : vertex_)
(*ME())(thel,fhel,ahel) +=
vert->evaluate(scale,wave_[ahel],
wavebar_[fhel],tensors_[thel]);
}
else {
(*ME())(thel,ahel,fhel) = 0.;
for(auto vert : vertex_)
(*ME())(thel,ahel,fhel) +=
vert->evaluate(scale,wave_[ahel],
wavebar_[fhel],tensors_[thel]);
}
}
}
}
double output = (ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
// colour and identical particle factors
- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
- decay[1]->dataPtr());
+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
// return the answer
return output;
}
Energy TFFDecayer::partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
if( inpart.second < outa.second + outb.second ) return ZERO;
if(perturbativeVertex_.size()==1 &&
perturbativeVertex_[0]) {
Energy2 scale = sqr(inpart.second);
tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
perturbativeVertex_[0]->setCoupling(scale, in, outa.first, outb.first);
double musq = sqr(outa.second/inpart.second);
double b = sqrt(1- 4.*musq);
double me2 = b*b*(5-2*b*b)*scale/120.*UnitRemoval::InvE2;
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
outb.second);
Energy output = norm(perturbativeVertex_[0]->norm())*me2*pcm/(8.*Constants::pi);
// colour factor
output *= colourFactor(inpart.first,outa.first,outb.first);
// return the answer
return output;
}
else {
- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
}
}
double TFFDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
// work out which is the fermion and antifermion
int ianti(0), iferm(1), iglu(2);
int itype[2];
for(unsigned int ix=0;ix<2;++ix) {
if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
else itype[ix] = 2;
}
if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(meopt==Initialize) {
// create tensor wavefunction for decaying particle
TensorWaveFunction::
calculateWaveFunctions(tensors3_, rho3_, const_ptr_cast<tPPtr>(&inpart), incoming, false);
}
// setup spin information when needed
if(meopt==Terminate) {
TensorWaveFunction::
constructSpinInfo(tensors3_, const_ptr_cast<tPPtr>(&inpart),incoming,true, false);
SpinorBarWaveFunction::
constructSpinInfo(wavebar3_ ,decay[iferm],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
VectorWaveFunction::
constructSpinInfo(gluon_ ,decay[iglu ],outgoing,true,false);
return 0.;
}
// calculate colour factors and number of colour flows
unsigned int nflow;
vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);
vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin2, PDT::Spin1Half,
PDT::Spin1Half, PDT::Spin1)));
// create wavefunctions
SpinorBarWaveFunction::
calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
SpinorWaveFunction::
calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
VectorWaveFunction::
calculateWaveFunctions(gluon_ , decay[iglu ],outgoing,true);
// gauge invariance test
#ifdef GAUGE_CHECK
gluon_.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) gluon_.push_back(VectorWaveFunction());
else {
gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
decay[iglu ]->dataPtr(),10,
outgoing));
}
}
#endif
if (! (outgoingVertex1_[inter] && outgoingVertex2_[inter]))
throw Exception()
<< "Invalid vertices for QCD radiation in TFF decay in TFFDecayer::threeBodyME"
<< Exception::runerror;
// identify fermion and/or anti-fermion vertex
AbstractFFVVertexPtr outgoingVertexF = outgoingVertex1_[inter];
AbstractFFVVertexPtr outgoingVertexA = outgoingVertex2_[inter];
if(outgoingVertex1_[inter]!=outgoingVertex2_[inter] &&
outgoingVertex1_[inter]->isIncoming(getParticleData(decay[ianti]->id())))
swap (outgoingVertexF, outgoingVertexA);
if(! (inpart.dataPtr()->iColour()==PDT::Colour0)){
throw Exception()
<< "Invalid vertices for QCD radiation in TFF decay in TFFDecayer::threeBodyME"
<< Exception::runerror;
}
Energy2 scale(sqr(inpart.mass()));
- const GeneralTwoBodyDecayer::CFlow & colourFlow
+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
= colourFlows(inpart, decay);
double gs(0.);
bool couplingSet(false);
#ifdef GAUGE_CHECK
double total=0.;
#endif
for(unsigned int it = 0; it < 5; ++it) {
for(unsigned int ifm = 0; ifm < 2; ++ifm) {
for(unsigned int ia = 0; ia < 2; ++ia) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from outgoing fermion
if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexF);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iferm]->dataPtr();
if(off->CC()) off = off->CC();
SpinorBarWaveFunction interS =
outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
gluon_[2*ig],decay[iferm]->mass());
assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ia], interS,tensors3_[it]);
if(!couplingSet) {
gs = abs(outgoingVertexF->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
(*ME[colourFlow[1][ix].first])(it, ifm, ia, ig) +=
colourFlow[1][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing antifermion
if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexA);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[ianti]->dataPtr();
if(off->CC()) off = off->CC();
SpinorWaveFunction interS =
outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
gluon_[2*ig],decay[ianti]->mass());
assert(wave3_[ia].particle()->id()==interS.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,interS,wavebar3_[ifm],tensors3_[it]);
if(!couplingSet) {
gs = abs(outgoingVertexA->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
(*ME[colourFlow[2][ix].first])(it, ifm, ia, ig) +=
colourFlow[2][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from 4 point vertex
if (fourPointVertex_[inter]) {
Complex diag = fourPointVertex_[inter]->evaluate(scale, wave3_[ia], wavebar3_[ifm],
gluon_[2*ig], tensors3_[it]);
for(unsigned int ix=0;ix<colourFlow[3].size();++ix) {
(*ME[colourFlow[3][ix].first])(it, ifm, ia, ig) +=
colourFlow[3][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
}
}
}
}
// contract matrices
double output=0.;
for(unsigned int ix=0; ix<nflow; ++ix){
for(unsigned int iy=0; iy<nflow; ++iy){
output+=cfactors[ix][iy]*(ME[ix]->contract(*ME[iy],rho3_)).real();
}
}
// divide by alpha_(s,em)
output *= (4.*Constants::pi)/sqr(gs);
#ifdef GAUGE_CHECK
double ratio = output/total;
if(abs(ratio)>1e-20) {
generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
for(unsigned int ix=0;ix<decay.size();++ix)
generator()->log() << *decay[ix] << "\n";
generator()->log() << "Test of gauge invariance " << ratio << "\n";
}
#endif
// return the answer
return output;
}
diff --git a/Decay/General/TFFDecayer.h b/Decay/General/TFFDecayer.h
--- a/Decay/General/TFFDecayer.h
+++ b/Decay/General/TFFDecayer.h
@@ -1,223 +1,232 @@
// -*- C++ -*-
//
// TFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_TFFDecayer_H
#define HERWIG_TFFDecayer_H
//
// This is the declaration of the TFFDecayer class.
//
-#include "GeneralTwoBodyDecayer.h"
+#include "GeneralTwoBodyDecayer2.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Helicity/Vertex/Tensor/FFTVertex.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Helicity/Vertex/Tensor/FFVTVertex.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::FFTVertexPtr;
/** \ingroup Decay
* The TFFDecayer class implements the decay of a tensor
* to 2 fermions in a general model. It holds an FFTVertex pointer
* that must be typecast from the VertexBase pointer held in
- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
* partialWidth().
*
- * @see GeneralTwoBodyDecayer
+ * @see GeneralTwoBodyDecayer2
*/
-class TFFDecayer: public GeneralTwoBodyDecayer {
+class TFFDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
TFFDecayer() {}
/** @name Virtual functions required by the Decayer class. */
//@{
/**
* 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 ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
- * @param decay The particles produced in the decay.
- * @param meopt Option for the matrix element
+ * @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.
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ 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;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa One of the decay products.
* @param outb The other decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {
POWHEGType output = FSR;
for(auto vertex : vertex_) {
if(vertex->orderInAllCouplings()!=1) {
output = No;
break;
}
}
return output;
}
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt);
/**
* Set the information on the decay
*/
virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr>,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > &,
map<ShowerInteraction,VertexBasePtr>);
//@}
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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
TFFDecayer & operator=(const TFFDecayer &);
private:
/**
* Abstract pointer to AbstractFFTVertex
*/
vector<AbstractFFTVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<FFTVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;
/**
* Abstract pointer to AbstractFFVTVertex for QCD radiation from 4 point vertex
*/
map<ShowerInteraction,AbstractFFVTVertexPtr> fourPointVertex_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Polarization tensors for the decaying particle
*/
mutable vector<TensorWaveFunction> tensors_;
/**
* Spinors for the decay products
*/
mutable vector<SpinorWaveFunction> wave_;
/**
* Barred spinors for the decay products
*/
mutable vector<SpinorBarWaveFunction> wavebar_;
/**
* Spin density matrix for 3 body decay
*/
mutable RhoDMatrix rho3_;
/**
* Tensor wavefunction for 3 body decay
*/
mutable vector<TensorWaveFunction> tensors3_;
/**
* Spinor wavefunction for 3 body decay
*/
mutable vector<SpinorWaveFunction> wave3_;
/**
* Barred spinor wavefunction for 3 body decay
*/
mutable vector<SpinorBarWaveFunction> wavebar3_;
/**
* Vector wavefunction for 3 body decay
*/
mutable vector<VectorWaveFunction> gluon_;
};
}
#endif /* HERWIG_TFFDecayer_H */
diff --git a/Decay/General/VFFDecayer.cc b/Decay/General/VFFDecayer.cc
--- a/Decay/General/VFFDecayer.cc
+++ b/Decay/General/VFFDecayer.cc
@@ -1,472 +1,475 @@
// -*- C++ -*-
//
// VFFDecayer.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VFFDecayer class.
//
#include "VFFDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Helicity/WaveFunction/VectorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorWaveFunction.h"
#include "ThePEG/Helicity/WaveFunction/SpinorBarWaveFunction.h"
#include "Herwig/Utilities/Kinematics.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
IBPtr VFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr VFFDecayer::fullclone() const {
return new_ptr(*this);
}
void VFFDecayer::setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr> vertex,
map<ShowerInteraction,VertexBasePtr> & inV,
const vector<map<ShowerInteraction,VertexBasePtr> > & outV,
map<ShowerInteraction,VertexBasePtr> ) {
decayInfo(incoming,outgoing);
for(auto vert : vertex) {
vertex_ .push_back(dynamic_ptr_cast<AbstractFFVVertexPtr>(vert));
perturbativeVertex_.push_back(dynamic_ptr_cast<FFVVertexPtr> (vert));
}
vector<ShowerInteraction> itemp={ShowerInteraction::QCD,ShowerInteraction::QED};
for(auto & inter : itemp) {
incomingVertex_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(inV.at(inter));
outgoingVertex1_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
outgoingVertex2_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
}
}
void VFFDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << perturbativeVertex_
<< incomingVertex_ << outgoingVertex1_
<< outgoingVertex2_;
}
void VFFDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> perturbativeVertex_
>> incomingVertex_ >> outgoingVertex1_
>> outgoingVertex2_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<VFFDecayer,GeneralTwoBodyDecayer>
+DescribeClass<VFFDecayer,GeneralTwoBodyDecayer2>
describeHerwigVFFDecayer("Herwig::VFFDecayer", "Herwig.so");
void VFFDecayer::Init() {
static ClassDocumentation<VFFDecayer> documentation
("The VFFDecayer implements the matrix element for the"
" decay of a vector to fermion-antifermion pair");
}
-double VFFDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,
+void VFFDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ int iferm(1),ianti(0);
+ if(decay[0]->id()>0) swap(iferm,ianti);
+ VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&part),
+ incoming,true,false);
+ SpinorBarWaveFunction::
+ constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
+ SpinorWaveFunction::
+ constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
+}
+
+double VFFDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
int iferm(1),ianti(0);
- if(decay[0]->id()>0) swap(iferm,ianti);
+ if(outgoing[0]->id()>0) swap(iferm,ianti);
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1,PDT::Spin1Half,PDT::Spin1Half)));
if(meopt==Initialize) {
VectorWaveFunction::calculateWaveFunctions(vectors_,rho_,
- const_ptr_cast<tPPtr>(&inpart),
- incoming,false);
+ const_ptr_cast<tPPtr>(&part),
+ incoming,false);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- VectorWaveFunction::constructSpinInfo(vectors_,const_ptr_cast<tPPtr>(&inpart),
- incoming,true,false);
- SpinorBarWaveFunction::
- constructSpinInfo(wavebar_,decay[iferm],outgoing,true);
- SpinorWaveFunction::
- constructSpinInfo(wave_ ,decay[ianti],outgoing,true);
- return 0.;
- }
SpinorBarWaveFunction::
- calculateWaveFunctions(wavebar_,decay[iferm],outgoing);
+ calculateWaveFunctions(wavebar_,momenta[iferm],outgoing[iferm],Helicity::outgoing);
SpinorWaveFunction::
- calculateWaveFunctions(wave_ ,decay[ianti],outgoing);
+ calculateWaveFunctions(wave_ ,momenta[ianti],outgoing[ianti],Helicity::outgoing);
// compute the matrix element
- Energy2 scale(inpart.mass()*inpart.mass());
+ Energy2 scale(part.mass()*part.mass());
for(unsigned int ifm = 0; ifm < 2; ++ifm) { //loop over fermion helicities
for(unsigned int ia = 0; ia < 2; ++ia) {// loop over antifermion helicities
for(unsigned int vhel = 0; vhel < 3; ++vhel) {//loop over vector helicities
- if(iferm > ianti) {
- (*ME())(vhel, ia, ifm) = 0.;
- for(auto vert : vertex_)
- (*ME())(vhel, ia, ifm) +=
- vert->evaluate(scale,wave_[ia],
- wavebar_[ifm],vectors_[vhel]);
- }
- else {
- (*ME())(vhel,ifm,ia)= 0.;
- for(auto vert : vertex_)
- (*ME())(vhel,ifm,ia) +=
- vert->evaluate(scale,wave_[ia],
- wavebar_[ifm],vectors_[vhel]);
- }
+ if(iferm > ianti) {
+ (*ME())(vhel, ia, ifm) = 0.;
+ for(auto vert : vertex_)
+ (*ME())(vhel, ia, ifm) +=
+ vert->evaluate(scale,wave_[ia],
+ wavebar_[ifm],vectors_[vhel]);
+ }
+ else {
+ (*ME())(vhel,ifm,ia)= 0.;
+ for(auto vert : vertex_)
+ (*ME())(vhel,ifm,ia) +=
+ vert->evaluate(scale,wave_[ia],
+ wavebar_[ifm],vectors_[vhel]);
+ }
}
}
}
double output=(ME()->contract(rho_)).real()/scale*UnitRemoval::E2;
// colour and identical particle factors
- output *= colourFactor(inpart.dataPtr(),decay[0]->dataPtr(),
- decay[1]->dataPtr());
+ output *= colourFactor(part.dataPtr(),outgoing[0],outgoing[1]);
// return the answer
return output;
}
Energy VFFDecayer::partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
if( inpart.second < outa.second + outb.second ) return ZERO;
if(perturbativeVertex_.size()==1 &&
perturbativeVertex_[0]) {
double mu1(outa.second/inpart.second), mu2(outb.second/inpart.second);
tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), outa.first, outb.first,in);
Complex cl(perturbativeVertex_[0]->left()), cr(perturbativeVertex_[0]->right());
double me2 = (norm(cl) + norm(cr))*( sqr(sqr(mu1) - sqr(mu2))
+ sqr(mu1) + sqr(mu2) - 2.)
- 6.*(cl*conj(cr) + cr*conj(cl)).real()*mu1*mu2;
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,outa.second,
outb.second);
Energy output = -norm(perturbativeVertex_[0]->norm())*me2*pcm /
(24.*Constants::pi);
// colour factor
output *= colourFactor(inpart.first,outa.first,outb.first);
// return the answer
return output;
}
else {
- return GeneralTwoBodyDecayer::partialWidth(inpart,outa,outb);
+ return GeneralTwoBodyDecayer2::partialWidth(inpart,outa,outb);
}
}
double VFFDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
// work out which is the fermion and antifermion
int ianti(0), iferm(1), iglu(2);
int itype[2];
for(unsigned int ix=0;ix<2;++ix) {
if(decay[ix]->dataPtr()->CC()) itype[ix] = decay[ix]->id()>0 ? 0:1;
else itype[ix] = 2;
}
if(itype[0]==0 && itype[1]!=0) swap(iferm, ianti);
if(itype[0]==2 && itype[1]==1) swap(iferm, ianti);
if(itype[0]==0 && itype[1]==0 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(itype[0]==1 && itype[1]==1 && decay[0]->dataPtr()->id()<decay[1]->dataPtr()->id())
swap(iferm, ianti);
if(meopt==Initialize) {
// create vector wavefunction for decaying particle
VectorWaveFunction::calculateWaveFunctions(vector3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
incoming, false);
}
// setup spin information when needed
if(meopt==Terminate) {
VectorWaveFunction::
constructSpinInfo(vector3_ ,const_ptr_cast<tPPtr>(&inpart),outgoing,true,false);
SpinorBarWaveFunction::
constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
SpinorWaveFunction::
constructSpinInfo(wave3_ ,decay[ianti],outgoing,true);
VectorWaveFunction::
constructSpinInfo(gluon_ ,decay[iglu ],outgoing,true,false);
return 0.;
}
// calculate colour factors and number of colour flows
unsigned int nflow;
vector<DVector> cfactors = getColourFactors(inpart, decay, nflow);
vector<GeneralDecayMEPtr> ME(nflow,new_ptr(GeneralDecayMatrixElement(PDT::Spin1, PDT::Spin1Half,
PDT::Spin1Half, PDT::Spin1)));
// create wavefunctions
SpinorBarWaveFunction::
calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
SpinorWaveFunction::
calculateWaveFunctions(wave3_ , decay[ianti],outgoing);
VectorWaveFunction::
calculateWaveFunctions(gluon_ , decay[iglu ],outgoing,true);
// gauge invariance test
#ifdef GAUGE_CHECK
gluon_.clear();
for(unsigned int ix=0;ix<3;++ix) {
if(ix==1) gluon_.push_back(VectorWaveFunction());
else {
gluon_.push_back(VectorWaveFunction(decay[iglu ]->momentum(),
decay[iglu ]->dataPtr(),10,
outgoing));
}
}
#endif
// identify fermion and/or anti-fermion vertex
AbstractFFVVertexPtr outgoingVertexF;
AbstractFFVVertexPtr outgoingVertexA;
identifyVertices(iferm, ianti, inpart, decay, outgoingVertexF, outgoingVertexA,inter);
Energy2 scale(sqr(inpart.mass()));
- const GeneralTwoBodyDecayer::CFlow & colourFlow
+ const GeneralTwoBodyDecayer2::CFlow & colourFlow
= colourFlows(inpart, decay);
double gs(0.);
bool couplingSet(false);
#ifdef GAUGE_CHECK
double total=0.;
#endif
for(unsigned int iv = 0; iv < 3; ++iv) {
for(unsigned int ifm = 0; ifm < 2; ++ifm) {
for(unsigned int ia = 0; ia < 2; ++ia) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from the incoming vector
if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(incomingVertex_[inter]);
VectorWaveFunction vectorInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),vector3_[iv],
gluon_[2*ig],inpart.mass());
assert(vector3_[iv].particle()->id()==vectorInter.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ia],wavebar3_[ifm],vectorInter);
if(!couplingSet) {
gs = abs(incomingVertex_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
(*ME[colourFlow[0][ix].first])(iv, ia, ifm, ig) +=
colourFlow[0][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing fermion
if((decay[iferm]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[iferm]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexF);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iferm]->dataPtr();
if(off->CC()) off = off->CC();
SpinorBarWaveFunction interS =
outgoingVertexF->evaluate(scale,3,off,wavebar3_[ifm],
gluon_[2*ig],decay[iferm]->mass());
assert(wavebar3_[ifm].particle()->id()==interS.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ia], interS,vector3_[iv]);
if(!couplingSet) {
gs = abs(outgoingVertexF->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[1].size();++ix) {
(*ME[colourFlow[1][ix].first])(iv, ia, ifm, ig) +=
colourFlow[1][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing antifermion
if((decay[ianti]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[ianti]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexA);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[ianti]->dataPtr();
if(off->CC()) off = off->CC();
SpinorWaveFunction interS =
outgoingVertexA->evaluate(scale,3,off,wave3_[ia],
gluon_[2*ig],decay[ianti]->mass());
assert(wave3_[ia].particle()->id()==interS.particle()->id());
Complex diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,interS,wavebar3_[ifm],vector3_[iv]);
if(!couplingSet) {
gs = abs(outgoingVertexA->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[2].size();++ix) {
(*ME[colourFlow[2][ix].first])(iv, ia, ifm, ig) +=
colourFlow[2][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
}
}
}
}
// contract matrices
double output=0.;
for(unsigned int ix=0; ix<nflow; ++ix){
for(unsigned int iy=0; iy<nflow; ++iy){
output+=cfactors[ix][iy]*(ME[ix]->contract(*ME[iy],rho3_)).real();
}
}
// divide by alpha_(S,EM)
output*=(4.*Constants::pi)/sqr(gs);
#ifdef GAUGE_CHECK
double ratio = output/total;
if(abs(ratio)>1e-20) {
generator()->log() << "Test of gauge invariance in decay\n" << inpart << "\n";
for(unsigned int ix=0;ix<decay.size();++ix)
generator()->log() << *decay[ix] << "\n";
generator()->log() << "Test of gauge invariance " << ratio << "\n";
}
#endif
//return output
return output;
}
void VFFDecayer::identifyVertices(const int iferm, const int ianti,
const Particle & inpart, const ParticleVector & decay,
AbstractFFVVertexPtr & outgoingVertexF,
AbstractFFVVertexPtr & outgoingVertexA,
ShowerInteraction inter){
// QCD vertices
if(inter==ShowerInteraction::QCD) {
// work out which fermion each outgoing vertex corresponds to
// two outgoing vertices
if( inpart.dataPtr() ->iColour()==PDT::Colour0 &&
((decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar) ||
(decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour8))){
if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
else if(inpart.dataPtr() ->iColour()==PDT::Colour8 &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
if(outgoingVertex1_[inter]==outgoingVertex2_[inter]){
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else if (outgoingVertex2_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
// one outgoing vertex
else if(inpart.dataPtr()->iColour()==PDT::Colour3){
if(decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour0){
if (outgoingVertex1_[inter]) outgoingVertexF = outgoingVertex1_[inter];
else if(outgoingVertex2_[inter]) outgoingVertexF = outgoingVertex2_[inter];
}
else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour8){
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr()))){
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
else {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
}
}
else if(inpart.dataPtr()->iColour()==PDT::Colour3bar){
if(decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour0){
if (outgoingVertex1_[inter]) outgoingVertexA = outgoingVertex1_[inter];
else if(outgoingVertex2_[inter]) outgoingVertexA = outgoingVertex2_[inter];
}
else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour8 &&
decay[ianti]->dataPtr()->iColour()==PDT::Colour3bar){
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))){
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
}
else if(inpart.dataPtr()->iColour()==PDT::Colour6 ||
inpart.dataPtr()->iColour()==PDT::Colour6bar) {
if (outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr()))) {
outgoingVertexF = outgoingVertex1_[inter];
outgoingVertexA = outgoingVertex2_[inter];
}
else {
outgoingVertexF = outgoingVertex2_[inter];
outgoingVertexA = outgoingVertex1_[inter];
}
}
if (! ((incomingVertex_[inter] && (outgoingVertexF || outgoingVertexA)) ||
( outgoingVertexF && outgoingVertexA)))
throw Exception()
<< "Invalid vertices for QCD radiation in VFF decay in VFFDecayer::identifyVertices"
<< Exception::runerror;
}
// QED
else {
if(decay[iferm]->dataPtr()->charged()) {
if (outgoingVertex1_[inter] &&
outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[iferm]->dataPtr())))
outgoingVertexF = outgoingVertex1_[inter];
else
outgoingVertexF = outgoingVertex2_[inter];
}
if(decay[ianti]->dataPtr()->charged()) {
if (outgoingVertex1_[inter] &&
outgoingVertex1_[inter]->isIncoming(const_ptr_cast<tPDPtr>(decay[ianti]->dataPtr())))
outgoingVertexA = outgoingVertex1_[inter];
else
outgoingVertexA = outgoingVertex2_[inter];
}
}
}
diff --git a/Decay/General/VFFDecayer.h b/Decay/General/VFFDecayer.h
--- a/Decay/General/VFFDecayer.h
+++ b/Decay/General/VFFDecayer.h
@@ -1,232 +1,241 @@
// -*- C++ -*-
//
// VFFDecayer.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_VFFDecayer_H
#define HERWIG_VFFDecayer_H
//
// This is the declaration of the VFFDecayer class.
//
-#include "GeneralTwoBodyDecayer.h"
+#include "GeneralTwoBodyDecayer2.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Helicity/Vertex/Vector/FFVVertex.h"
#include "ThePEG/Helicity/Vertex/Vector/VVVVertex.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::FFVVertexPtr;
/** \ingroup Decay
* The VFFDecayer class implements the decay of a vector
* to 2 fermions in a general model. It holds an FFVVertex pointer
* that must be typecast from the VertexBase pointer held in
- * GeneralTwoBodyDecayer. It implents the virtual functions me2() and
+ * GeneralTwoBodyDecayer2. It implents the virtual functions me2() and
* partialWidth().
*
- * @see GeneralTwoBodyDecayer
+ * @see GeneralTwoBodyDecayer2
*/
-class VFFDecayer: public GeneralTwoBodyDecayer {
+class VFFDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
VFFDecayer() {}
public:
/** @name Virtual functions required by the Decayer class. */
//@{
/**
- * Return the matrix element squared for a given mode and phase-space channel.
- * @param ichan The channel we are calculating the matrix element for.
+ * Return the matrix element squared for a given mode and phase-space channel.
+ * @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
- * @param decay The particles produced in the decay.
- * @param meopt Option for the matrix element
+ * @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.
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ 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;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa One of the decay products.
* @param outb The other decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const;
/**
* Has a POWHEG style correction
*/
virtual POWHEGType hasPOWHEGCorrection() {
POWHEGType output = FSR;
for(auto vertex : vertex_) {
if(vertex->orderInAllCouplings()!=1) {
output = No;
break;
}
}
return output;
}
/**
* Three-body matrix element including additional QCD radiation
*/
virtual double threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt);
/**
* Indentify outgoing vertices for the fermion and antifermion
*/
void identifyVertices(const int iferm, const int ianti,
const Particle & inpart, const ParticleVector & decay,
AbstractFFVVertexPtr & abstractOutgoingVertexF,
AbstractFFVVertexPtr & abstractOutgoingVertexA,
ShowerInteraction inter);
/**
* Set the information on the decay
*/
virtual void setDecayInfo(PDPtr incoming, PDPair outgoing,
vector<VertexBasePtr>,
map<ShowerInteraction,VertexBasePtr> &,
const vector<map<ShowerInteraction,VertexBasePtr> > &,
map<ShowerInteraction,VertexBasePtr>);
//@}
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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
VFFDecayer & operator=(const VFFDecayer &);
private:
/**
* Abstract pointer to AbstractFFVVertex
*/
vector<AbstractFFVVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<FFVVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractVVVVertex for QCD radiation from incoming vector
*/
map<ShowerInteraction,AbstractVVVVertexPtr> incomingVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex1_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertex2_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Polarization vectors for the decaying particle
*/
mutable vector<VectorWaveFunction> vectors_;
/**
* Spinors for the decay products
*/
mutable vector<SpinorWaveFunction> wave_;
/**
* Barred spinors for the decay products
*/
mutable vector<SpinorBarWaveFunction> wavebar_;
/**
* Spin density matrix for 3 body decay
*/
mutable RhoDMatrix rho3_;
/**
* Scalar wavefunction for 3 body decay
*/
mutable vector<VectorWaveFunction> vector3_;
/**
* Spinor wavefunction for 3 body decay
*/
mutable vector<SpinorWaveFunction> wave3_;
/**
* Barred spinor wavefunction for 3 body decay
*/
mutable vector<SpinorBarWaveFunction> wavebar3_;
/**
* Vector wavefunction for 3 body decay
*/
mutable vector<VectorWaveFunction> gluon_;
};
}
#endif /* HERWIG_VFFDecayer_H */

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 8:16 PM (8 h, 57 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4226349
Default Alt Text
(95 KB)

Event Timeline