Page MenuHomeHEPForge
Files F8308999

No OneTemporary
Actions

View Options

diff --git a/Decay/General/FFSDecayer.cc b/Decay/General/FFSDecayer.cc
--- a/Decay/General/FFSDecayer.cc
+++ b/Decay/General/FFSDecayer.cc
@@ -1,465 +1,472 @@
// -*- C++ -*-
//
// FFSDecayer.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 FFSDecayer class.
//
#include "FFSDecayer.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 FFSDecayer::clone() const {
return new_ptr(*this);
}
IBPtr FFSDecayer::fullclone() const {
return new_ptr(*this);
}
void FFSDecayer::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<AbstractFFVVertexPtr>(inV.at(inter));
outgoingVertexF_[inter] = AbstractFFVVertexPtr();
outgoingVertexS_[inter] = AbstractVSSVertexPtr();
if(outV[0].at(inter)) {
if (outV[0].at(inter)->getName()==VertexType::FFV)
outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[0].at(inter));
else
outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[0].at(inter));
}
if(outV[1].at(inter)) {
if (outV[1].at(inter)->getName()==VertexType::FFV)
outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
else
outgoingVertexS_[inter] = dynamic_ptr_cast<AbstractVSSVertexPtr>(outV[1].at(inter));
}
}
}
void FFSDecayer::persistentOutput(PersistentOStream & os) const {
os << perturbativeVertex_ << vertex_
<< incomingVertex_ << outgoingVertexF_
<< outgoingVertexS_;
}
void FFSDecayer::persistentInput(PersistentIStream & is, int) {
is >> perturbativeVertex_ >> vertex_
>> incomingVertex_ >> outgoingVertexF_
>> outgoingVertexS_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<FFSDecayer,GeneralTwoBodyDecayer>
+DescribeClass<FFSDecayer,GeneralTwoBodyDecayer2>
describeHerwigFFSDecayer("Herwig::FFSDecayer", "Herwig.so");
void FFSDecayer::Init() {
static ClassDocumentation<FFSDecayer> documentation
("The FFSDecayer class implements the decay of a fermion to "
"a fermion and a scalar.");
}
-double FFSDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,
+void FFSDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ int itype[2];
+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0? 0:1;
+ else itype[0] = 2;
+ if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0? 0:1;
+ else itype[1] = 2;
+ bool ferm(itype[0] == 0 || itype[1] == 0 || (itype[0] == 2 && itype[1] == 2));
+ // for the decaying particle
+ if(ferm) {
+ SpinorWaveFunction::
+ constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&part),incoming,true);
+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
+ }
+ else {
+ SpinorBarWaveFunction::
+ constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&part),incoming,true);
+ SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
+ }
+ ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
+}
+
+double FFSDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin0)));
//Need to use different barred or unbarred spinors depending on
//whether particle is cc or not.
int itype[2];
- if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0? 0:1;
- else itype[0] = 2;
- if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0? 0:1;
+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0? 0:1;
+ else itype[0] = 2;
+ if(outgoing[0]->CC()) itype[1] = outgoing[0]->id() > 0? 0:1;
else itype[1] = 2;
bool ferm(itype[0] == 0 || itype[1] == 0 || (itype[0] == 2 && itype[1] == 2));
if(meopt==Initialize) {
// spinors and rho
if(ferm) {
SpinorWaveFunction ::calculateWaveFunctions(wave_,rho_,
- const_ptr_cast<tPPtr>(&inpart),
+ const_ptr_cast<tPPtr>(&part),
incoming);
if(wave_[0].wave().Type() != SpinorType::u)
for(unsigned int ix = 0; ix < 2; ++ix) wave_ [ix].conjugate();
}
else {
SpinorBarWaveFunction::calculateWaveFunctions(wavebar_,rho_,
- const_ptr_cast<tPPtr>(&inpart),
+ const_ptr_cast<tPPtr>(&part),
incoming);
if(wavebar_[0].wave().Type() != SpinorType::v)
for(unsigned int ix = 0; ix < 2; ++ix) wavebar_[ix].conjugate();
}
// fix rho if no correlations
fixRho(rho_);
}
- // setup spin info when needed
- if(meopt==Terminate) {
- // for the decaying particle
- if(ferm) {
- SpinorWaveFunction::
- constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
- SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
- }
- else {
- SpinorBarWaveFunction::
- constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
- SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
- }
- ScalarWaveFunction::constructSpinInfo(decay[1],outgoing,true);
- }
if(ferm)
SpinorBarWaveFunction::
- calculateWaveFunctions(wavebar_,decay[0],outgoing);
+ calculateWaveFunctions(wavebar_,momenta[0],outgoing[0],Helicity::outgoing);
else
SpinorWaveFunction::
- calculateWaveFunctions(wave_ ,decay[0],outgoing);
- ScalarWaveFunction scal(decay[1]->momentum(),decay[1]->dataPtr(),outgoing);
- Energy2 scale(sqr(inpart.mass()));
+ calculateWaveFunctions(wave_ ,momenta[0],outgoing[0],Helicity::outgoing);
+ ScalarWaveFunction scal(momenta[1],outgoing[1],Helicity::outgoing);
+ Energy2 scale(sqr(part.mass()));
for(unsigned int if1 = 0; if1 < 2; ++if1) {
for(unsigned int if2 = 0; if2 < 2; ++if2) {
if(ferm) (*ME())(if1, if2, 0) = 0.;
else (*ME())(if2, if1, 0) = 0.;
for(auto vert : vertex_) {
if(ferm) (*ME())(if1, if2, 0) +=
vert->evaluate(scale,wave_[if1],wavebar_[if2],scal);
else (*ME())(if2, if1, 0) +=
vert->evaluate(scale,wave_[if1],wavebar_[if2],scal);
}
}
}
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 FFSDecayer::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(0.),mu2(0.);
tcPDPtr in = inpart.first->CC() ? tcPDPtr(inpart.first->CC()) : inpart.first;
if(outa.first->iSpin() == PDT::Spin1Half) {
mu1 = outa.second/inpart.second;
mu2 = outb.second/inpart.second;
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in, outa.first, outb.first);
}
else {
mu1 = outb.second/inpart.second;
mu2 = outa.second/inpart.second;
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in, outb.first, outa.first);
}
double c2 = norm(perturbativeVertex_[0]->norm());
Complex cl = perturbativeVertex_[0]->left();
Complex cr = perturbativeVertex_[0]->right();
double me2 = c2*( (norm(cl) + norm(cr))*(1. + sqr(mu1) - sqr(mu2))
+ 2.*mu1*(conj(cl)*cr + conj(cr)*cl).real() );
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
outb.second);
Energy output = me2*pcm/16./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 FFSDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
int iscal (0), iferm (1), iglu (2);
// get location of outgoing fermion/scalar
if(decay[1]->dataPtr()->iSpin()==PDT::Spin0) swap(iscal,iferm);
// work out whether inpart is a fermion or antifermion
int itype[2];
if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
else itype[0] = 2;
if(decay[iferm]->dataPtr()->CC()) itype[1] = decay[iferm]->id() > 0 ? 0 : 1;
else itype[1] = 2;
bool ferm(false);
if(itype[0] == itype[1] ) {
ferm = itype[0]==0 || (itype[0]==2 && decay[iscal]->id() < 0);
}
else if(itype[0] == 2) {
ferm = itype[1]==0;
}
else if(itype[1] == 2) {
ferm = itype[0]==0;
}
else if((itype[0] == 1 && itype[1] == 0) ||
(itype[0] == 0 && itype[1] == 1)) {
if(abs(inpart.id())<=16) {
ferm = itype[0]==0;
}
else if(abs(decay[iferm]->id())<=16) {
ferm = itype[1]==0;
}
else {
ferm = true;
}
}
else
assert(false);
if(meopt==Initialize) {
// create spinor (bar) for decaying particle
if(ferm) {
SpinorWaveFunction::calculateWaveFunctions(wave3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
incoming);
if(wave3_[0].wave().Type() != SpinorType::u)
for(unsigned int ix = 0; ix < 2; ++ix) wave3_[ix].conjugate();
}
else {
SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_,rho3_, const_ptr_cast<tPPtr>(&inpart),
incoming);
if(wavebar3_[0].wave().Type() != SpinorType::v)
for(unsigned int ix = 0; ix < 2; ++ix) wavebar3_[ix].conjugate();
}
}
// setup spin information when needed
if(meopt==Terminate) {
if(ferm) {
SpinorWaveFunction::
constructSpinInfo(wave3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorBarWaveFunction::constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
}
else {
SpinorBarWaveFunction::
constructSpinInfo(wavebar3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorWaveFunction::constructSpinInfo(wave3_,decay[iferm],outgoing,true);
}
ScalarWaveFunction::constructSpinInfo( decay[iscal],outgoing,true);
VectorWaveFunction::constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
return 0.;
}
// calulate 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::Spin1Half, PDT::Spin0,
PDT::Spin1Half, PDT::Spin1)));
// create wavefunctions
if (ferm) SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
else SpinorWaveFunction:: calculateWaveFunctions(wave3_ , decay[iferm],outgoing);
ScalarWaveFunction swave3_(decay[iscal]->momentum(), decay[iscal]->dataPtr(),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 (! ((incomingVertex_[inter] && (outgoingVertexF_[inter] || outgoingVertexS_[inter])) ||
(outgoingVertexF_[inter] && outgoingVertexS_[inter])))
throw Exception()
<< "Invalid vertices for radiation in FFS decay in FFSDecayer::threeBodyME"
<< Exception::runerror;
// sort out colour flows
int F(1), S(2);
if (decay[iscal]->dataPtr()->iColour()==PDT::Colour3 &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour8)
swap(F,S);
else if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[iscal]->dataPtr()->iColour()==PDT::Colour8)
swap(F,S);
Complex diag;
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 ifi = 0; ifi < 2; ++ifi) {
for(unsigned int ifo = 0; ifo < 2; ++ifo) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from the incoming fermion
if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(incomingVertex_[inter]);
if (ferm) {
SpinorWaveFunction spinorInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wave3_[ifi],
gluon_[2*ig],inpart.mass());
assert(wave3_[ifi].particle()->id()==spinorInter.particle()->id());
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifo],swave3_);
}
else {
SpinorBarWaveFunction spinorBarInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wavebar3_[ifi],
gluon_[2*ig],inpart.mass());
assert(wavebar3_[ifi].particle()->id()==spinorBarInter.particle()->id());
diag = 0.;
for(auto vertex :vertex_)
diag+= vertex->evaluate(scale,wave3_[ifo], spinorBarInter,swave3_);
}
if(!couplingSet) {
gs = abs(incomingVertex_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
(*ME[colourFlow[0][ix].first])(ifi, 0, ifo, 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_[inter]);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iferm]->dataPtr();
if(off->CC()) off = off->CC();
if (ferm) {
SpinorBarWaveFunction spinorBarInter =
outgoingVertexF_[inter]->evaluate(scale,3,off,wavebar3_[ifo],
gluon_[2*ig],decay[iferm]->mass());
assert(wavebar3_[ifo].particle()->id()==spinorBarInter.particle()->id());
diag = 0.;
for(auto vertex :vertex_)
diag+= vertex->evaluate(scale,wave3_[ifi],spinorBarInter,swave3_);
}
else {
SpinorWaveFunction spinorInter =
outgoingVertexF_[inter]->evaluate(scale,3,off,wave3_[ifo],
gluon_[2*ig],decay[iferm]->mass());
assert(wave3_[ifo].particle()->id()==spinorInter.particle()->id());
diag = 0.;
for(auto vertex :vertex_)
diag+= vertex->evaluate(scale,spinorInter,wavebar3_[ifi],swave3_);
}
if(!couplingSet) {
gs = abs(outgoingVertexF_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[F].size();++ix) {
(*ME[colourFlow[F][ix].first])(ifi, 0, ifo, ig) +=
colourFlow[F][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
// radiation from outgoing scalar
if((decay[iscal]->dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(decay[iscal]->dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(outgoingVertexS_[inter]);
// ensure you get correct ougoing 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],
swave3_,decay[iscal]->mass());
assert(swave3_.particle()->id()==scalarInter.particle()->id());
if (ferm){
diag = 0.;
for(auto vertex :vertex_)
diag += vertex->evaluate(scale,wave3_[ifi],wavebar3_[ifo],scalarInter);
}
else {
diag = 0.;
for(auto vertex :vertex_)
diag += vertex->evaluate(scale,wave3_[ifo],wavebar3_[ifi],scalarInter);
}
if(!couplingSet) {
gs = abs(outgoingVertexS_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[S].size();++ix) {
(*ME[colourFlow[S][ix].first])(ifi, 0, ifo, ig) +=
colourFlow[S][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/FFSDecayer.h b/Decay/General/FFSDecayer.h
--- a/Decay/General/FFSDecayer.h
+++ b/Decay/General/FFSDecayer.h
@@ -1,218 +1,227 @@
// -*- C++ -*-
//
// FFSDecayer.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_FFSDecayer_H
#define HERWIG_FFSDecayer_H
//
// This is the declaration of the FFSDecayer 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 FFSDecayer class implements the decay of a fermion
* to a fermion and a vector 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 FFSDecayer: public GeneralTwoBodyDecayer {
+class FFSDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
FFSDecayer() {}
/** @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);
/**
* 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.
*/
FFSDecayer & operator=(const FFSDecayer &);
private:
/**
* Abstract pointer to AbstractFFSVertex
*/
vector<AbstractFFSVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<FFSVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from incoming (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> incomingVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertexF_;
/**
* Abstract pointer to AbstractVSSVertex for QCD radiation from outgoing scalar
*/
map<ShowerInteraction,AbstractVSSVertexPtr> outgoingVertexS_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Spinor wavefunctions
*/
mutable vector<SpinorWaveFunction> wave_ ;
/**
* Barred spinor wavefunctions
*/
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_FFSDecayer_H */
diff --git a/Decay/General/FFVDecayer.cc b/Decay/General/FFVDecayer.cc
--- a/Decay/General/FFVDecayer.cc
+++ b/Decay/General/FFVDecayer.cc
@@ -1,458 +1,469 @@
// -*- C++ -*-
//
// FFVDecayer.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 FFVDecayer class.
//
#include "FFVDecayer.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 "Herwig/Decay/GeneralDecayMatrixElement.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
IBPtr FFVDecayer::clone() const {
return new_ptr(*this);
}
IBPtr FFVDecayer::fullclone() const {
return new_ptr(*this);
}
void FFVDecayer::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<AbstractFFVVertexPtr>(inV.at(inter));
if(outV[0].at(inter)) {
if (outV[0].at(inter)->getName()==VertexType::FFV)
outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(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::FFV)
outgoingVertexF_[inter] = dynamic_ptr_cast<AbstractFFVVertexPtr>(outV[1].at(inter));
else
outgoingVertexV_[inter] = dynamic_ptr_cast<AbstractVVVVertexPtr>(outV[1].at(inter));
}
}
}
void FFVDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << perturbativeVertex_
<< incomingVertex_ << outgoingVertexF_
<< outgoingVertexV_;
}
void FFVDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> perturbativeVertex_
>> incomingVertex_ >> outgoingVertexF_
>> outgoingVertexV_;
}
-double FFVDecayer::me2(const int , const Particle & inpart,
- const ParticleVector & decay,
+void FFVDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ // type of process
+ int itype[2];
+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0 ? 0 : 1;
+ else itype[0] = 2;
+ if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0 ? 0 : 1;
+ else itype[1] = 2;
+ //Need to use different barred or unbarred spinors depending on
+ //whether particle is cc or not.
+ bool ferm(itype[0] == 0 || itype[1] == 0 || (itype[0] == 2 && itype[1] == 2));
+ // for the decaying particle
+ if(ferm) {
+ SpinorWaveFunction::
+ constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&part),incoming,true);
+ SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
+ }
+ else {
+ SpinorBarWaveFunction::
+ constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&part),incoming,true);
+ SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
+ }
+ VectorWaveFunction::
+ constructSpinInfo(vector_,decay[1],outgoing,true,false);
+}
+
+double FFVDecayer::me2(const int,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
if(!ME())
ME(new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,PDT::Spin1)));
// type of process
int itype[2];
- if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
- else itype[0] = 2;
- if(decay[0]->dataPtr()->CC()) itype[1] = decay[0]->id() > 0 ? 0 : 1;
- else itype[1] = 2;
+ if(part.dataPtr()->CC()) itype[0] = part.id() > 0 ? 0 : 1;
+ else itype[0] = 2;
+ if(outgoing[0]->CC()) itype[1] = outgoing[0]->id() > 0 ? 0 : 1;
+ else itype[1] = 2;
//Need to use different barred or unbarred spinors depending on
//whether particle is cc or not.
bool ferm(itype[0] == 0 || itype[1] == 0 || (itype[0] == 2 && itype[1] == 2));
if(meopt==Initialize) {
// spinors and rho
if(ferm) {
SpinorWaveFunction ::calculateWaveFunctions(wave_,rho_,
- const_ptr_cast<tPPtr>(&inpart),
+ const_ptr_cast<tPPtr>(&part),
incoming);
if(wave_[0].wave().Type() != SpinorType::u)
for(unsigned int ix = 0; ix < 2; ++ix) wave_ [ix].conjugate();
}
else {
SpinorBarWaveFunction::calculateWaveFunctions(wavebar_,rho_,
- const_ptr_cast<tPPtr>(&inpart),
+ const_ptr_cast<tPPtr>(&part),
incoming);
if(wavebar_[0].wave().Type() != SpinorType::v)
for(unsigned int ix = 0; ix < 2; ++ix) wavebar_[ix].conjugate();
}
// fix rho if no correlations
fixRho(rho_);
}
- // setup spin info when needed
- if(meopt==Terminate) {
- // for the decaying particle
- if(ferm) {
- SpinorWaveFunction::
- constructSpinInfo(wave_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
- SpinorBarWaveFunction::constructSpinInfo(wavebar_,decay[0],outgoing,true);
- }
- else {
- SpinorBarWaveFunction::
- constructSpinInfo(wavebar_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
- SpinorWaveFunction::constructSpinInfo(wave_,decay[0],outgoing,true);
- }
- VectorWaveFunction::
- constructSpinInfo(vector_,decay[1],outgoing,true,false);
- }
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
if(ferm)
SpinorBarWaveFunction::
- calculateWaveFunctions(wavebar_,decay[0],outgoing);
+ calculateWaveFunctions(wavebar_,momenta[0],outgoing[0],Helicity::outgoing);
else
SpinorWaveFunction::
- calculateWaveFunctions(wave_ ,decay[0],outgoing);
- bool massless = decay[1]->dataPtr()->mass()==ZERO;
+ calculateWaveFunctions(wave_ ,momenta[0],outgoing[0],Helicity::outgoing);
+ bool massless = outgoing[1]->mass()==ZERO;
VectorWaveFunction::
- calculateWaveFunctions(vector_,decay[1],outgoing,massless);
+ calculateWaveFunctions(vector_,momenta[1],outgoing[1],Helicity::outgoing,massless);
for(unsigned int if1 = 0; if1 < 2; ++if1) {
for(unsigned int if2 = 0; if2 < 2; ++if2) {
for(unsigned int vhel = 0; vhel < 3; ++vhel) {
if(massless && vhel == 1) ++vhel;
if(ferm)
(*ME())(if1, if2,vhel) = 0.;
else
(*ME())(if2, if1, vhel) = 0.;
for(auto vertex : vertex_) {
if(ferm)
(*ME())(if1, if2,vhel) +=
vertex->evaluate(scale,wave_[if1],wavebar_[if2],vector_[vhel]);
else
(*ME())(if2, if1, vhel) +=
vertex->evaluate(scale,wave_[if1],wavebar_[if2],vector_[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 FFVDecayer::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;
if( outa.first->iSpin() == PDT::Spin1Half)
perturbativeVertex_[0]->setCoupling(sqr(inpart.second), in,
outa.first, outb.first);
else {
swap(mu1,mu2);
perturbativeVertex_[0]->setCoupling(sqr(inpart.second),in,
outb.first,outa.first);
}
Complex cl(perturbativeVertex_[0]->left()),cr(perturbativeVertex_[0]->right());
double me2(0.);
if( mu2 > 0. ) {
me2 = (norm(cl) + norm(cr))*(1. + sqr(mu1*mu2) + sqr(mu2)
- 2.*sqr(mu1) - 2.*sqr(mu2*mu2)
+ sqr(mu1*mu1))
- 6.*mu1*sqr(mu2)*(conj(cl)*cr + conj(cr)*cl).real();
me2 /= sqr(mu2);
}
else
me2 = 2.*( (norm(cl) + norm(cr))*(sqr(mu1) + 1.)
- 4.*mu1*(conj(cl)*cr + conj(cr)*cl).real() );
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second, outa.second,
outb.second);
Energy output = norm(perturbativeVertex_[0]->norm())*me2*pcm/16./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);
}
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeClass<FFVDecayer,GeneralTwoBodyDecayer>
+DescribeClass<FFVDecayer,GeneralTwoBodyDecayer2>
describeHerwigFFVDecayer("Herwig::FFVDecayer", "Herwig.so");
void FFVDecayer::Init() {
static ClassDocumentation<FFVDecayer> documentation
("The FFVDecayer class implements the decay of a fermion to a fermion and a vector boson");
}
double FFVDecayer::threeBodyME(const int , const Particle & inpart,
const ParticleVector & decay,
ShowerInteraction inter, MEOption meopt) {
int iferm (0), ivect (1), iglu (2);
// get location of outgoing lepton/vector
if(decay[1]->dataPtr()->iSpin()==PDT::Spin1Half) swap(iferm,ivect);
// work out whether inpart is a fermion or antifermion
int itype[2];
if(inpart.dataPtr()->CC()) itype[0] = inpart.id() > 0 ? 0 : 1;
else itype[0] = 2;
if(decay[iferm]->dataPtr()->CC()) itype[1] = decay[iferm]->id() > 0 ? 0 : 1;
else itype[1] = 2;
bool ferm(itype[0] == 0 || itype[1] == 0 ||
(itype[0] == 2 && itype[1] == 2 && decay[ivect]->id() < 0));
// no emissions from massive vectors
bool massless = decay[ivect]->dataPtr()->mass()==ZERO;
if(meopt==Initialize) {
// create spinor (bar) for decaying particle
if(ferm) {
SpinorWaveFunction::calculateWaveFunctions(wave3_, rho3_, const_ptr_cast<tPPtr>(&inpart),
incoming);
if(wave3_[0].wave().Type() != SpinorType::u)
for(unsigned int ix = 0; ix < 2; ++ix) wave3_[ix].conjugate();
}
else {
SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_,rho3_, const_ptr_cast<tPPtr>(&inpart),
incoming);
if(wavebar3_[0].wave().Type() != SpinorType::v)
for(unsigned int ix = 0; ix < 2; ++ix) wavebar3_[ix].conjugate();
}
}
// setup spin information when needed
if(meopt==Terminate) {
if(ferm) {
SpinorWaveFunction::
constructSpinInfo(wave3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorBarWaveFunction::constructSpinInfo(wavebar3_,decay[iferm],outgoing,true);
}
else {
SpinorBarWaveFunction::
constructSpinInfo(wavebar3_,const_ptr_cast<tPPtr>(&inpart),incoming,true);
SpinorWaveFunction::constructSpinInfo(wave3_,decay[iferm],outgoing,true);
}
VectorWaveFunction::constructSpinInfo(vector3_, decay[ivect],outgoing,true,massless);
VectorWaveFunction::constructSpinInfo(gluon_, decay[iglu ],outgoing,true,false);
return 0.;
}
// calulate 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::Spin1Half, PDT::Spin1Half,
PDT::Spin1, PDT::Spin1)));
// create wavefunctions
if (ferm) SpinorBarWaveFunction::calculateWaveFunctions(wavebar3_, decay[iferm],outgoing);
else SpinorWaveFunction:: calculateWaveFunctions(wave3_ , decay[iferm],outgoing);
VectorWaveFunction::calculateWaveFunctions(vector3_, decay[ivect],outgoing,massless);
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] && (outgoingVertexF_[inter] || outgoingVertexV_[inter])) ||
(outgoingVertexF_[inter] && outgoingVertexV_[inter])))
throw Exception()
<< "Invalid vertices for QCD radiation in FFV decay in FFVDecayer::threeBodyME"
<< Exception::runerror;
// sort out colour flows
int F(1), V(2);
if (decay[iferm]->dataPtr()->iColour()==PDT::Colour3bar &&
decay[ivect]->dataPtr()->iColour()==PDT::Colour8)
swap(F,V);
else if (decay[ivect]->dataPtr()->iColour()==PDT::Colour3 &&
decay[iferm]->dataPtr()->iColour()==PDT::Colour8)
swap(F,V);
Complex diag;
Energy2 scale(sqr(inpart.mass()));
- const GeneralTwoBodyDecayer::CFlow & colourFlow = colourFlows(inpart, decay);
+ const GeneralTwoBodyDecayer2::CFlow & colourFlow = colourFlows(inpart, decay);
double gs(0.);
bool couplingSet(false);
#ifdef GAUGE_CHECK
double total=0.;
#endif
for(unsigned int ifi = 0; ifi < 2; ++ifi) {
for(unsigned int ifo = 0; ifo < 2; ++ifo) {
for(unsigned int iv = 0; iv < 3; ++iv) {
for(unsigned int ig = 0; ig < 2; ++ig) {
// radiation from the incoming fermion
if((inpart.dataPtr()->coloured() && inter==ShowerInteraction::QCD) ||
(inpart.dataPtr()->charged() && inter==ShowerInteraction::QED) ) {
assert(incomingVertex_[inter]);
if (ferm){
SpinorWaveFunction spinorInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wave3_[ifi],
gluon_[2*ig],inpart.mass());
assert(wave3_[ifi].particle()->id()==spinorInter.particle()->id());
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifo],vector3_[iv]);
}
else {
SpinorBarWaveFunction spinorBarInter =
incomingVertex_[inter]->evaluate(scale,3,inpart.dataPtr(),wavebar3_[ifi],
gluon_[2*ig],inpart.mass());
assert(wavebar3_[ifi].particle()->id()==spinorBarInter.particle()->id());
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ifo], spinorBarInter,vector3_[iv]);
}
if(!couplingSet) {
gs = abs(incomingVertex_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[0].size();++ix) {
(*ME[colourFlow[0][ix].first])(ifi, ifo, iv, 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_[inter]);
// ensure you get correct outgoing particle from first vertex
tcPDPtr off = decay[iferm]->dataPtr();
if(off->CC()) off = off->CC();
if (ferm) {
SpinorBarWaveFunction spinorBarInter =
outgoingVertexF_[inter]->evaluate(scale,3,off,wavebar3_[ifo],
gluon_[2*ig],decay[iferm]->mass());
assert(wavebar3_[ifo].particle()->id()==spinorBarInter.particle()->id());
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ifi],spinorBarInter,vector3_[iv]);
}
else {
SpinorWaveFunction spinorInter =
outgoingVertexF_[inter]->evaluate(scale,3,off,wave3_[ifo],
gluon_[2*ig],decay[iferm]->mass());
assert(wave3_[ifo].particle()->id()==spinorInter.particle()->id());
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,spinorInter,wavebar3_[ifi],vector3_[iv]);
}
if(!couplingSet) {
gs = abs(outgoingVertexF_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[F].size();++ix) {
(*ME[colourFlow[F][ix].first])(ifi, ifo, iv, ig) +=
colourFlow[F][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 ougoing 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 (ferm) {
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ifi],wavebar3_[ifo],vectorInter);
}
else {
diag = 0.;
for(auto vertex : vertex_)
diag += vertex->evaluate(scale,wave3_[ifo],wavebar3_[ifi],vectorInter);
}
if(!couplingSet) {
gs = abs(outgoingVertexV_[inter]->norm());
couplingSet = true;
}
for(unsigned int ix=0;ix<colourFlow[V].size();++ix) {
(*ME[colourFlow[V][ix].first])(ifi, ifo, iv, ig) +=
colourFlow[V][ix].second*diag;
}
#ifdef GAUGE_CHECK
total+=norm(diag);
#endif
}
}
if(massless) ++iv;
}
}
}
// 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/FFVDecayer.h b/Decay/General/FFVDecayer.h
--- a/Decay/General/FFVDecayer.h
+++ b/Decay/General/FFVDecayer.h
@@ -1,223 +1,232 @@
// -*- C++ -*-
//
// FFVDecayer.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_FFVDecayer_H
#define HERWIG_FFVDecayer_H
//
// This is the declaration of the FFVDecayer 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 FFVDecayer class implements the decay of a fermion
* to a fermion and a vector 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 FFVDecayer: public GeneralTwoBodyDecayer {
+class FFVDecayer: public GeneralTwoBodyDecayer2 {
public:
/**
* The default constructor.
*/
FFVDecayer() {}
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);
/**
* 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.
*/
FFVDecayer & operator=(const FFVDecayer &);
private:
/**
* Abstract pointer to AbstractFFVVertex
*/
vector<AbstractFFVVertexPtr> vertex_;
/**
* Pointer to the perturbative vertex
*/
vector<FFVVertexPtr> perturbativeVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from incoming (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> incomingVertex_;
/**
* Abstract pointer to AbstractFFVVertex for QCD radiation from outgoing (anti)fermion
*/
map<ShowerInteraction,AbstractFFVVertexPtr> outgoingVertexF_;
/**
* Abstract pointer to AbstractVVVVertex for QCD radiation from outgoing vector
*/
map<ShowerInteraction,AbstractVVVVertexPtr> outgoingVertexV_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Spinor wavefunction
*/
mutable vector<SpinorWaveFunction> wave_ ;
/**
* Barred spinor wavefunction
*/
mutable vector<SpinorBarWaveFunction> wavebar_;
/**
* Polarization vectors
*/
mutable vector<VectorWaveFunction> vector_;
/**
* Spin density matrix for 3 body decay
*/
mutable RhoDMatrix rho3_;
/**
* vector 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 gluon in 3 body decay
*/
mutable vector<VectorWaveFunction> gluon_;
};
}
#endif /* HERWIG_FFVDecayer_H */

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 1:54 PM (11 h, 52 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4018828
Default Alt Text
(50 KB)

Event Timeline