Page MenuHomeHEPForge
Files F7877610

No OneTemporary
Actions

View Options

diff --git a/Decay/General/FtoFFFDecayer.cc b/Decay/General/FtoFFFDecayer.cc
--- a/Decay/General/FtoFFFDecayer.cc
+++ b/Decay/General/FtoFFFDecayer.cc
@@ -1,311 +1,313 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FtoFFFDecayer class.
//
#include "FtoFFFDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
IBPtr FtoFFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr FtoFFFDecayer::fullclone() const {
return new_ptr(*this);
}
void FtoFFFDecayer::persistentOutput(PersistentOStream & os) const {
os << sca_ << vec_ << ten_;
}
void FtoFFFDecayer::persistentInput(PersistentIStream & is, int) {
is >> sca_ >> vec_ >> ten_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<FtoFFFDecayer,GeneralThreeBodyDecayer>
describeHerwigFtoFFFDecayer("Herwig::FtoFFFDecayer", "Herwig.so");
void FtoFFFDecayer::Init() {
static ClassDocumentation<FtoFFFDecayer> documentation
("The FtoFFFDecayer class implements the general decay of a fermion to "
"three fermions.");
}
void FtoFFFDecayer::doinit() {
GeneralThreeBodyDecayer::doinit();
if(outgoing().empty()) return;
unsigned int ndiags = getProcessInfo().size();
sca_.resize(ndiags);
vec_.resize(ndiags);
ten_.resize(ndiags);
for(unsigned int ix = 0;ix < ndiags; ++ix) {
TBDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if( offshell->CC() ) offshell = offshell->CC();
if(offshell->iSpin() == PDT::Spin0) {
AbstractFFSVertexPtr vert1 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.first);
AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in FtoFFFDecayer::doinit()"
<< Exception::runerror;
sca_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractFFVVertexPtr vert1 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a vector diagram in FtoFFFDecayer::doinit()"
<< Exception::runerror;
vec_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin2) {
AbstractFFTVertexPtr vert1 = dynamic_ptr_cast<AbstractFFTVertexPtr>
(current.vertices.first);
AbstractFFTVertexPtr vert2 = dynamic_ptr_cast<AbstractFFTVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a tensor diagram in FtoFFFDecayer::doinit()"
<< Exception::runerror;
ten_[ix] = make_pair(vert1, vert2);
}
}
}
-double FtoFFFDecayer::me2(const int ichan, const Particle & inpart,
- const ParticleVector & decay,
+void FtoFFFDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ // for the decaying particle
+ if(part.id()<0)
+ SpinorWaveFunction::constructSpinInfo(inwave_.first,
+ const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ else
+ SpinorBarWaveFunction::constructSpinInfo(inwave_.second,
+ const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ // outgoing particles
+ for(unsigned int ix = 0; ix < 3; ++ix) {
+ SpinorWaveFunction::
+ constructSpinInfo(outwave_[ix].first,decay[ix],Helicity::outgoing,true);
+ }
+}
+
+double FtoFFFDecayer::me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
// particle or CC of particle
- bool cc = (*getProcessInfo().begin()).incoming != inpart.id();
+ bool cc = (*getProcessInfo().begin()).incoming != part.id();
// special handling or first/last call
const vector<vector<double> > cfactors(getColourFactors());
const vector<vector<double> > nfactors(getLargeNcColourFactors());
const size_t ncf(numberOfFlows());
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
if(meopt==Initialize) {
SpinorWaveFunction::
- calculateWaveFunctions(inwave_.first,rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(inwave_.first,rho_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming);
inwave_.second.resize(2);
if(inwave_.first[0].wave().Type() == SpinorType::u) {
for(unsigned int ix = 0; ix < 2; ++ix) {
inwave_.second[ix] = inwave_.first[ix].bar();
inwave_.second[ix].conjugate();
}
}
else {
for(unsigned int ix = 0; ix < 2; ++ix) {
inwave_.second[ix] = inwave_.first[ix].bar();
inwave_.first[ix].conjugate();
}
}
// fix rho if no correlations
fixRho(rho_);
}
- // setup spin info when needed
- if(meopt==Terminate) {
- // for the decaying particle
- if(inpart.id()<0)
- SpinorWaveFunction::constructSpinInfo(inwave_.first,
- const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- else
- SpinorBarWaveFunction::constructSpinInfo(inwave_.second,
- const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- // outgoing particles
- for(unsigned int ix = 0; ix < 3; ++ix) {
- SpinorWaveFunction::
- constructSpinInfo(outwave_[ix].first,decay[ix],Helicity::outgoing,true);
- }
- }
// outgoing particles
for(unsigned int ix = 0; ix < 3; ++ix) {
SpinorWaveFunction::
- calculateWaveFunctions(outwave_[ix].first,decay[ix],Helicity::outgoing);
+ calculateWaveFunctions(outwave_[ix].first,momenta[ix],outgoing[ix],Helicity::outgoing);
outwave_[ix].second.resize(2);
if(outwave_[ix].first[0].wave().Type() == SpinorType::u) {
for(unsigned int iy = 0; iy < 2; ++iy) {
outwave_[ix].second[iy] = outwave_[ix].first[iy].bar();
outwave_[ix].first[iy].conjugate();
}
}
else {
for(unsigned int iy = 0; iy < 2; ++iy) {
outwave_[ix].second[iy] = outwave_[ix].first[iy].bar();
outwave_[ix].second[iy].conjugate();
}
}
}
- bool ferm = inpart.id()>0;
+ bool ferm = part.id()>0;
vector<Complex> flows(ncf, Complex(0.)),largeflows(ncf, Complex(0.));
static const unsigned int out2[3]={1,0,0},out3[3]={2,2,1};
vector<GeneralDecayMEPtr> mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half)));
vector<GeneralDecayMEPtr> mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,PDT::Spin1Half,
PDT::Spin1Half,PDT::Spin1Half)));
unsigned int ihel[4];
for(ihel[0] = 0; ihel[0] < 2; ++ihel[0]) {
for(ihel[1] = 0; ihel[1] < 2; ++ihel[1]) {
for(ihel[2] = 0; ihel[2] < 2; ++ihel[2]) {
for(ihel[3] = 0; ihel[3] < 2; ++ihel[3]) {
flows = vector<Complex>(ncf, Complex(0.));
largeflows = vector<Complex>(ncf, Complex(0.));
unsigned int idiag=0;
for(vector<TBDiagram>::const_iterator dit=getProcessInfo().begin();
dit!=getProcessInfo().end();++dit) {
if(ichan>=0&&diagramMap()[ichan]!=idiag) {
++idiag;
continue;
}
// the sign from normal ordering
double sign = ferm ? 1. : -1;
// outgoing wavefunction and NO sign
if (dit->channelType==TBDiagram::channel23) sign *= -1.;
else if(dit->channelType==TBDiagram::channel13) sign *= 1.;
else if(dit->channelType==TBDiagram::channel12) sign *= -1.;
else throw Exception()
<< "Unknown diagram type in FtoFFFDecayer::me2()" << Exception::runerror;
// wavefunctions
SpinorWaveFunction w0,w3;
SpinorBarWaveFunction w1,w2;
// incoming wavefunction
if(ferm) {
w0 = inwave_.first [ihel[0]];
w1 = outwave_[dit->channelType].second[ihel[dit->channelType+1]];
}
else {
w0 = outwave_[dit->channelType].first [ihel[dit->channelType+1]];
w1 = inwave_.second[ihel[0]];
}
- if(decay[out2[dit->channelType]]->id()<0&&
- decay[out3[dit->channelType]]->id()>0) {
+ if(outgoing[out2[dit->channelType]]->id()<0&&
+ outgoing[out3[dit->channelType]]->id()>0) {
w2 = outwave_[out3[dit->channelType]].second[ihel[out3[dit->channelType]+1]];
w3 = outwave_[out2[dit->channelType]].first [ihel[out2[dit->channelType]+1]];
sign *= -1.;
}
else {
w2 = outwave_[out2[dit->channelType]].second[ihel[out2[dit->channelType]+1]];
w3 = outwave_[out3[dit->channelType]].first [ihel[out3[dit->channelType]+1]];
}
tcPDPtr offshell = dit->intermediate;
if(cc&&offshell->CC()) offshell=offshell->CC();
Complex diag(0.);
// intermediate scalar
if (offshell->iSpin() == PDT::Spin0) {
ScalarWaveFunction inters = sca_[idiag].first->
evaluate(scale, widthOption(), offshell, w0, w1);
diag = sca_[idiag].second->evaluate(scale,w3,w2,inters);
}
// intermediate vector
else if(offshell->iSpin() == PDT::Spin1) {
VectorWaveFunction interv = vec_[idiag].first->
evaluate(scale, widthOption(), offshell, w0, w1);
diag = vec_[idiag].second->evaluate(scale,w3,w2,interv);
}
// intermediate tensor
else if(offshell->iSpin() == PDT::Spin2) {
TensorWaveFunction intert = ten_[idiag].first->
evaluate(scale, widthOption(), offshell, w0, w1);
diag = ten_[idiag].second->evaluate(scale,w3,w2,intert);
}
// unknown
else throw Exception()
<< "Unknown intermediate in FtoFFFDecayer::me2()"
<< Exception::runerror;
// apply NO sign
diag *= sign;
// matrix element for the different colour flows
if(ichan<0) {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
else {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
++idiag;
}
// now add the flows to the me2 with appropriate colour factors
for(unsigned int ix = 0; ix < ncf; ++ix) {
(*mes[ix])(ihel[0],ihel[1],ihel[2],ihel[3]) = flows[ix];
(*mel[ix])(ihel[0],ihel[1],ihel[2],ihel[3]) = largeflows[ix];
}
}
}
}
}
double me2(0.);
if(ichan<0) {
vector<double> pflows(ncf,0.);
for(unsigned int ix = 0; ix < ncf; ++ix) {
for(unsigned int iy = 0; iy < ncf; ++ iy) {
double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real();
me2 += con;
if(ix==iy) {
con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real();
pflows[ix] += con;
}
}
}
double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.));
ptotal *=UseRandom::rnd();
for(unsigned int ix=0;ix<pflows.size();++ix) {
if(ptotal<=pflows[ix]) {
colourFlow(ix);
ME(mes[ix]);
break;
}
ptotal-=pflows[ix];
}
}
else {
unsigned int iflow = colourFlow();
me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real();
}
// return the matrix element squared
return me2;
}
WidthCalculatorBasePtr FtoFFFDecayer::
threeBodyMEIntegrator(const DecayMode & ) const {
vector<int> intype;
vector<Energy> inmass,inwidth;
vector<double> inpow,inweights;
constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights);
return new_ptr(ThreeBodyAllOnCalculator<FtoFFFDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,0,
outgoing()[0]->mass(),outgoing()[1]->mass(),outgoing()[2]->mass(),
relativeError()));
}
diff --git a/Decay/General/FtoFFFDecayer.h b/Decay/General/FtoFFFDecayer.h
--- a/Decay/General/FtoFFFDecayer.h
+++ b/Decay/General/FtoFFFDecayer.h
@@ -1,142 +1,151 @@
// -*- C++ -*-
#ifndef HERWIG_FtoFFFDecayer_H
#define HERWIG_FtoFFFDecayer_H
//
// This is the declaration of the FtoFFFDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* The FtoFFFDecayer class provides the general matrix elements for the
* decay of a (anti)fermion to three (anti)fermions.
*
* @see \ref FtoFFFDecayerInterfaces "The interfaces"
* defined for FtoFFFDecayer.
*/
class FtoFFFDecayer: public GeneralThreeBodyDecayer {
public:
+
+ /**
+ * Return the matrix element squared for a given mode and phase-space channel.
+ * @param ichan The channel we are calculating the matrix element for.
+ * @param part The decaying Particle.
+ * @param outgoing The particles produced in the decay
+ * @param momenta The momenta of the particles produced in the decay
+ * @param meopt Option for the calculation of the matrix element
+ * @return The matrix element squared for the phase-space configuration.
+ */
+ double me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ MEOption meopt) const;
/**
- * 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
- * @return The matrix element squared for the phase-space configuration.
+ * Construct the SpinInfos for the particles produced in the decay
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ virtual void constructSpinInfo(const Particle & part,
+ ParticleVector outgoing) const;
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/ void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FtoFFFDecayer & operator=(const FtoFFFDecayer &);
private:
/**
* Store the vector of FFSVertex pairs
*/
vector<pair<AbstractFFSVertexPtr, AbstractFFSVertexPtr> > sca_;
/**
* Store the vector of FFVVertex pairs
*/
vector<pair<AbstractFFVVertexPtr, AbstractFFVVertexPtr> > vec_;
/**
* Store the vector of FFTVertex pairs
*/
vector<pair<AbstractFFTVertexPtr, AbstractFFTVertexPtr> > ten_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Spinors for incoming particle
*/
mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > inwave_;
/**
* Spinors for outgoing particles
*/
mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > outwave_[3];
};
}
#endif /* HERWIG_FtoFFFDecayer_H */
diff --git a/Decay/General/FtoFVVDecayer.cc b/Decay/General/FtoFVVDecayer.cc
--- a/Decay/General/FtoFVVDecayer.cc
+++ b/Decay/General/FtoFVVDecayer.cc
@@ -1,403 +1,403 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FtoFVVDecayer class.
//
#include "FtoFVVDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Decay/DecayPhaseSpaceMode.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
IBPtr FtoFVVDecayer::clone() const {
return new_ptr(*this);
}
IBPtr FtoFVVDecayer::fullclone() const {
return new_ptr(*this);
}
void FtoFVVDecayer::persistentOutput(PersistentOStream & os) const {
os << sca_ << fer_ << vec_ << ten_;
}
void FtoFVVDecayer::persistentInput(PersistentIStream & is, int) {
is >> sca_ >> fer_ >> vec_ >> ten_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<FtoFVVDecayer,GeneralThreeBodyDecayer>
describeHerwigFtoFVVDecayer("Herwig::FtoFVVDecayer", "Herwig.so");
void FtoFVVDecayer::Init() {
static ClassDocumentation<FtoFVVDecayer> documentation
("The FtoFVVDecayer class implements the general decay of a fermion to "
"a fermion and a pair of vectors.");
}
void FtoFVVDecayer::doinit() {
GeneralThreeBodyDecayer::doinit();
if(outgoing().empty()) return;
unsigned int ndiags = getProcessInfo().size();
sca_.resize(ndiags);
fer_.resize(ndiags);
vec_.resize(ndiags);
ten_.resize(ndiags);
for(unsigned int ix = 0;ix < ndiags; ++ix) {
TBDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if(offshell->iSpin() == PDT::Spin0) {
AbstractFFSVertexPtr vert1 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.first);
AbstractVVSVertexPtr vert2 = dynamic_ptr_cast<AbstractVVSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in FtoFVVDecayer::doinit()"
<< Exception::runerror;
sca_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1Half) {
AbstractFFVVertexPtr vert1 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in FtoFVVDecayer::doinit()"
<< Exception::runerror;
fer_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractFFVVertexPtr vert1 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.first);
AbstractVVVVertexPtr vert2 = dynamic_ptr_cast<AbstractVVVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a vector diagram in FtoFVVDecayer::doinit()"
<< Exception::runerror;
vec_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin2) {
AbstractFFTVertexPtr vert1 = dynamic_ptr_cast<AbstractFFTVertexPtr>
(current.vertices.first);
AbstractVVTVertexPtr vert2 = dynamic_ptr_cast<AbstractVVTVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a tensor diagram in FtoFVVDecayer::doinit()"
<< Exception::runerror;
ten_[ix] = make_pair(vert1, vert2);
}
}
}
+void FtoFVVDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ bool ferm( part.id() > 0 );
+ // for the decaying particle
+ if(ferm)
+ SpinorWaveFunction::constructSpinInfo(fwave_,
+ const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ else
+ SpinorBarWaveFunction::constructSpinInfo(fbwave_,
+ const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ int ivec(-1);
+ // outgoing particles
+ for(int ix = 0; ix < 3; ++ix) {
+ tPPtr p = decay[ix];
+ if( p->dataPtr()->iSpin() == PDT::Spin1Half ) {
+ if( ferm ) {
+ SpinorBarWaveFunction::
+ constructSpinInfo(fbwave_, p, Helicity::outgoing,true);
+ }
+ else {
+ SpinorWaveFunction::
+ constructSpinInfo(fwave_ , p, Helicity::outgoing,true);
+ }
+ }
+ else if( p->dataPtr()->iSpin() == PDT::Spin1 ) {
+ if( ivec < 0 ) {
+ ivec = ix;
+ VectorWaveFunction::
+ constructSpinInfo(vwave_.first , p, Helicity::outgoing, true, false);
+ }
+ else {
+ VectorWaveFunction::
+ constructSpinInfo(vwave_.second, p, Helicity::outgoing, true, false);
+ }
+ }
+ }
+}
-double FtoFVVDecayer::me2(const int ichan, const Particle & inpart,
- const ParticleVector & decay,
+double FtoFVVDecayer::me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
// particle or CC of particle
- bool cc = (*getProcessInfo().begin()).incoming != inpart.id();
+ bool cc = (*getProcessInfo().begin()).incoming != part.id();
// special handling or first/last call
//Set up wave-functions
- bool ferm( inpart.id() > 0 );
+ bool ferm( part.id() > 0 );
if(meopt==Initialize) {
if( ferm ) {
SpinorWaveFunction::
- calculateWaveFunctions(fwave_,rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(fwave_,rho_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming);
if( fwave_[0].wave().Type() != SpinorType::u )
fwave_[0].conjugate();
if( fwave_[1].wave().Type() != SpinorType::u )
fwave_[1].conjugate();
}
else {
SpinorBarWaveFunction::
- calculateWaveFunctions(fbwave_, rho_, const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(fbwave_, rho_, const_ptr_cast<tPPtr>(&part),
Helicity::incoming);
if( fbwave_[0].wave().Type() != SpinorType::v )
fbwave_[0].conjugate();
if( fbwave_[1].wave().Type() != SpinorType::v )
fbwave_[1].conjugate();
}
// fix rho if no correlations
fixRho(rho_);
}
- // setup spin info when needed
- if(meopt==Terminate) {
- // for the decaying particle
- if(ferm)
- SpinorWaveFunction::constructSpinInfo(fwave_,
- const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- else
- SpinorBarWaveFunction::constructSpinInfo(fbwave_,
- const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- int ivec(-1);
- // outgoing particles
- for(int ix = 0; ix < 3; ++ix) {
- tPPtr p = decay[ix];
- if( p->dataPtr()->iSpin() == PDT::Spin1Half ) {
- if( ferm ) {
- SpinorBarWaveFunction::
- constructSpinInfo(fbwave_, p, Helicity::outgoing,true);
- }
- else {
- SpinorWaveFunction::
- constructSpinInfo(fwave_ , p, Helicity::outgoing,true);
- }
- }
- else if( p->dataPtr()->iSpin() == PDT::Spin1 ) {
- if( ivec < 0 ) {
- ivec = ix;
- VectorWaveFunction::
- constructSpinInfo(vwave_.first , p, Helicity::outgoing, true, false);
- }
- else {
- VectorWaveFunction::
- constructSpinInfo(vwave_.second, p, Helicity::outgoing, true, false);
- }
- }
- }
- return 0.;
- }
// outgoing, keep track of fermion and first occurrence of vector positions
int isp(-1), ivec(-1);
// outgoing particles
pair<bool,bool> mass = make_pair(false,false);
for(int ix = 0; ix < 3; ++ix) {
- tPPtr p = decay[ix];
- if( p->dataPtr()->iSpin() == PDT::Spin1Half ) {
+ if( outgoing[ix]->iSpin() == PDT::Spin1Half ) {
isp = ix;
if( ferm ) {
SpinorBarWaveFunction::
- calculateWaveFunctions(fbwave_, p, Helicity::outgoing);
+ calculateWaveFunctions(fbwave_, momenta[ix],outgoing[ix], Helicity::outgoing);
if( fbwave_[0].wave().Type() != SpinorType::u )
fbwave_[0].conjugate();
if( fbwave_[1].wave().Type() != SpinorType::u )
fbwave_[1].conjugate();
}
else {
SpinorWaveFunction::
- calculateWaveFunctions(fwave_, p, Helicity::outgoing);
+ calculateWaveFunctions(fwave_, momenta[ix],outgoing[ix], Helicity::outgoing);
if( fwave_[0].wave().Type() != SpinorType::v )
fwave_[0].conjugate();
if( fwave_[1].wave().Type() != SpinorType::v )
fwave_[1].conjugate();
}
}
- else if( p->dataPtr()->iSpin() == PDT::Spin1 ) {
- bool massless = p->id() == ParticleID::gamma || p->id() == ParticleID::g;
+ else if( outgoing[ix]->iSpin() == PDT::Spin1 ) {
+ bool massless = outgoing[ix]->id() == ParticleID::gamma || outgoing[ix]->id() == ParticleID::g;
if( ivec < 0 ) {
ivec = ix;
VectorWaveFunction::
- calculateWaveFunctions(vwave_.first , p, Helicity::outgoing, massless);
+ calculateWaveFunctions(vwave_.first , momenta[ix],outgoing[ix], Helicity::outgoing, massless);
mass.first = massless;
}
else {
VectorWaveFunction::
- calculateWaveFunctions(vwave_.second, p, Helicity::outgoing, massless);
+ calculateWaveFunctions(vwave_.second, momenta[ix],outgoing[ix], Helicity::outgoing, massless);
mass.second = massless;
}
}
}
assert(isp >= 0 && ivec >= 0);
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
const vector<vector<double> > cfactors(getColourFactors());
const vector<vector<double> > nfactors(getLargeNcColourFactors());
const size_t ncf(numberOfFlows());
vector<Complex> flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.));
vector<GeneralDecayMEPtr>
mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,
(isp == 0) ? PDT::Spin1Half : PDT::Spin1,
(isp == 1) ? PDT::Spin1Half : PDT::Spin1,
(isp == 2) ? PDT::Spin1Half : PDT::Spin1)));
vector<GeneralDecayMEPtr>
mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1Half,
(isp == 0) ? PDT::Spin1Half : PDT::Spin1,
(isp == 1) ? PDT::Spin1Half : PDT::Spin1,
(isp == 2) ? PDT::Spin1Half : PDT::Spin1)));
//Helicity calculation
for( unsigned int if1 = 0; if1 < 2; ++if1 ) {
for( unsigned int if2 = 0; if2 < 2; ++if2 ) {
for( unsigned int iv1 = 0; iv1 < 3; ++iv1 ) {
if ( mass.first && iv1 == 1 ) continue;
for( unsigned int iv2 = 0; iv2 < 3; ++iv2 ) {
if ( mass.second && iv2 == 1 ) continue;
flows = vector<Complex>(ncf, Complex(0.));
largeflows = vector<Complex>(ncf, Complex(0.));
unsigned int idiag(0);
for(vector<TBDiagram>::const_iterator dit = getProcessInfo().begin();
dit != getProcessInfo().end(); ++dit) {
// If we are selecting a particular channel
if( ichan >= 0 && diagramMap()[ichan] != idiag ) {
++idiag;
continue;
}
tcPDPtr offshell = (*dit).intermediate;
if(cc&&offshell->CC()) offshell=offshell->CC();
Complex diag;
if( offshell->iSpin() == PDT::Spin1Half ) {
// Make sure we connect the correct particles
VectorWaveFunction vw1, vw2;
if( (*dit).channelType == TBDiagram::channel23 ) {
vw1 = vwave_.first[iv1];
vw2 = vwave_.second[iv2];
}
else if( (*dit).channelType == TBDiagram::channel12 ) {
vw1 = vwave_.second[iv2];
vw2 = vwave_.first[iv1];
}
else {
if( ivec < isp ) {
vw1 = vwave_.second[iv2];
vw2 = vwave_.first[iv1];
}
else {
vw1 = vwave_.first[iv1];
vw2 = vwave_.second[iv2];
}
}
if( ferm ) {
SpinorWaveFunction inters =
fer_[idiag].first->evaluate(scale, widthOption(), offshell,
fwave_[if1], vw1);
diag = fer_[idiag].second->evaluate(scale, inters, fbwave_[if2],
vw2);
}
else {
SpinorBarWaveFunction inters =
fer_[idiag].first->evaluate(scale, widthOption(), offshell,
fbwave_[if2], vw1);
diag = fer_[idiag].second->evaluate(scale, fwave_[if1], inters,
vw2);
}
}
else if( offshell->iSpin() == PDT::Spin0 ) {
ScalarWaveFunction inters =
sca_[idiag].first->evaluate(scale, widthOption(), offshell,
fwave_[if1], fbwave_[if2]);
diag = sca_[idiag].second->evaluate(scale, vwave_.first[iv1],
vwave_.second[iv2], inters);
}
else if( offshell->iSpin() == PDT::Spin1 ) {
VectorWaveFunction interv =
vec_[idiag].first->evaluate(scale, widthOption(), offshell,
fwave_[if1], fbwave_[if2]);
diag = vec_[idiag].second->evaluate(scale, vwave_.first[iv1],
vwave_.second[iv2], interv);
}
else if( offshell->iSpin() == PDT::Spin2 ) {
TensorWaveFunction intert =
ten_[idiag].first->evaluate(scale, widthOption(), offshell,
fwave_[if1], fbwave_[if2]);
diag = ten_[idiag].second->evaluate(scale, vwave_.first[iv1],
vwave_.second[iv2], intert);
}
else
throw Exception()
<< "Unknown intermediate in FtoFVVDecayer::me2()"
<< Exception::runerror;
//NO sign
if( !ferm ) diag *= -1;
if(ichan<0) {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
else {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
++idiag;
}// end diagram loop
// now add the flows to the me2
unsigned int h1(if1), h2(if2);
if( !ferm ) swap(h1,h2);
for(unsigned int ix = 0; ix < ncf; ++ix) {
if(isp == 0) {
(*mes[ix])(h1, h2, iv1, iv2) = flows[ix];
(*mel[ix])(h1, h2, iv1, iv2) = largeflows[ix];
}
else if(isp == 1) {
(*mes[ix])(h1, iv1, h2, iv2) = flows[ix];
(*mel[ix])(h1, iv1, h2, iv2) = largeflows[ix];
}
else if(isp == 2) {
(*mes[ix])(h1, iv1, iv2, h2) = flows[ix];
(*mel[ix])(h1, iv1, h2, iv2) = largeflows[ix];
}
}
}//v2
}//v1
}//f2
}//f1
//Finally, work out me2. This depends on whether we are selecting channels
//or not
double me2(0.);
if(ichan<0) {
vector<double> pflows(ncf,0.);
for(unsigned int ix = 0; ix < ncf; ++ix) {
for(unsigned int iy = 0; iy < ncf; ++ iy) {
double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real();
me2 += con;
if(ix==iy) {
con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real();
pflows[ix] += con;
}
}
}
double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.));
ptotal *= UseRandom::rnd();
for(unsigned int ix=0;ix<pflows.size();++ix) {
if(ptotal<=pflows[ix]) {
colourFlow(ix);
ME(mes[ix]);
break;
}
ptotal-=pflows[ix];
}
}
else {
unsigned int iflow = colourFlow();
me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real();
}
// return the matrix element squared
return me2;
}
WidthCalculatorBasePtr FtoFVVDecayer::
threeBodyMEIntegrator(const DecayMode & ) const {
vector<int> intype;
vector<Energy> inmass,inwidth;
vector<double> inpow,inweights;
constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights);
return new_ptr(ThreeBodyAllOnCalculator<FtoFVVDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,0,
outgoing()[0]->mass(),outgoing()[1]->mass(),
outgoing()[2]->mass(),relativeError()));
}
diff --git a/Decay/General/FtoFVVDecayer.h b/Decay/General/FtoFVVDecayer.h
--- a/Decay/General/FtoFVVDecayer.h
+++ b/Decay/General/FtoFVVDecayer.h
@@ -1,156 +1,165 @@
// -*- C++ -*-
#ifndef HERWIG_FtoFVVDecayer_H
#define HERWIG_FtoFVVDecayer_H
//
// This is the declaration of the FtoFVVDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* The FtoFVVDecayer class provides the general matrix elements for the
* decay of a fermion to a fermion and two vector bosons.
*
* @see \ref FtoFVVDecayerInterfaces "The interfaces"
* defined for FtoFVVDecayer.
*/
class FtoFVVDecayer: public GeneralThreeBodyDecayer {
public:
-
+
/**
- * 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;
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FtoFVVDecayer & operator=(const FtoFVVDecayer &);
private:
/**
* Store the vector of scalar intermediates
*/
vector<pair<AbstractFFSVertexPtr, AbstractVVSVertexPtr> > sca_;
/**
* Store the vector for fermion intermediates
*/
vector<pair<AbstractFFVVertexPtr, AbstractFFVVertexPtr> > fer_;
/**
* Store the vector for gauge boson intermediates
*/
vector<pair<AbstractFFVVertexPtr, AbstractVVVVertexPtr> > vec_;
/**
* Store the vector of tensor intermediates
*/
vector<pair<AbstractFFTVertexPtr, AbstractVVTVertexPtr> > ten_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Spinor wavefunctions
*/
mutable vector<SpinorWaveFunction> fwave_;
/**
* Barred spinor wavefunctions
*/
mutable vector<SpinorBarWaveFunction> fbwave_;
/**
* Vector wavefunctions
*/
mutable pair<vector<VectorWaveFunction>, vector<VectorWaveFunction> > vwave_;
};
}
#endif /* HERWIG_FtoFVVDecayer_H */
diff --git a/Decay/General/GeneralThreeBodyDecayer.cc b/Decay/General/GeneralThreeBodyDecayer.cc
--- a/Decay/General/GeneralThreeBodyDecayer.cc
+++ b/Decay/General/GeneralThreeBodyDecayer.cc
@@ -1,1028 +1,1016 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GeneralThreeBodyDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
-#include "Herwig/Decay/DecayPhaseSpaceMode.h"
+#include "Herwig/Decay/PhaseSpaceMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
using namespace Herwig;
void GeneralThreeBodyDecayer::persistentOutput(PersistentOStream & os) const {
os << incoming_ << outgoing_ << diagrams_ << diagmap_ << colour_ << colourLargeNC_
<< nflow_ << widthOpt_ << refTag_ << refTagCC_ << intOpt_ << relerr_;
}
void GeneralThreeBodyDecayer::persistentInput(PersistentIStream & is, int) {
is >> incoming_ >> outgoing_ >> diagrams_ >> diagmap_ >> colour_ >> colourLargeNC_
>> nflow_ >> widthOpt_ >> refTag_ >> refTagCC_ >> intOpt_ >> relerr_;
}
// The following static variable is needed for the type
// description system in ThePEG.
-DescribeAbstractClass<GeneralThreeBodyDecayer,DecayIntegrator>
+DescribeAbstractClass<GeneralThreeBodyDecayer,DecayIntegrator2>
describeHerwigGeneralThreeBodyDecayer("Herwig::GeneralThreeBodyDecayer", "Herwig.so");
void GeneralThreeBodyDecayer::Init() {
static ClassDocumentation<GeneralThreeBodyDecayer> documentation
("The GeneralThreeBodyDecayer class is the base class for the implementation of"
" all three body decays based on spin structures in Herwig.");
static Switch<GeneralThreeBodyDecayer,unsigned int> interfaceWidthOption
("WidthOption",
"Option for the treatment of the widths of the intermediates",
&GeneralThreeBodyDecayer::widthOpt_, 1, false, false);
static SwitchOption interfaceWidthOptionFixed
(interfaceWidthOption,
"Fixed",
"Use fixed widths",
1);
static SwitchOption interfaceWidthOptionRunning
(interfaceWidthOption,
"Running",
"Use running widths",
2);
static SwitchOption interfaceWidthOptionZero
(interfaceWidthOption,
"Zero",
"Set the widths to zero",
3);
static Switch<GeneralThreeBodyDecayer,unsigned int> interfacePartialWidthIntegration
("PartialWidthIntegration",
"Switch to control the partial width integration",
&GeneralThreeBodyDecayer::intOpt_, 0, false, false);
static SwitchOption interfacePartialWidthIntegrationAllPoles
(interfacePartialWidthIntegration,
"AllPoles",
"Include all potential poles",
0);
static SwitchOption interfacePartialWidthIntegrationShallowestPole
(interfacePartialWidthIntegration,
"ShallowestPole",
"Only include the shallowest pole",
1);
static Parameter<GeneralThreeBodyDecayer,double> interfaceRelativeError
("RelativeError",
"The relative error for the GQ integration of the partial width",
&GeneralThreeBodyDecayer::relerr_, 1e-2, 1e-10, 1.,
false, false, Interface::limited);
}
ParticleVector GeneralThreeBodyDecayer::decay(const Particle & parent,
const tPDVector & children) const {
// return empty vector if products heavier than parent
Energy mout(ZERO);
for(tPDVector::const_iterator it=children.begin();
it!=children.end();++it) mout+=(**it).massMin();
if(mout>parent.mass()) return ParticleVector();
// generate the decay
bool cc;
int imode=modeNumber(cc,parent.dataPtr(),children);
// generate the kinematics
ParticleVector decay=generate(generateIntermediates(),cc,imode,parent);
// make the colour connections
colourConnections(parent, decay);
// return the answer
return decay;
}
int GeneralThreeBodyDecayer::
modeNumber(bool & cc, tcPDPtr in, const tPDVector & outin) const {
assert( !refTag_.empty() && !refTagCC_.empty() );
// check number of outgoing particles
if( outin.size() != 3 || abs(in->id()) != abs(incoming_->id()) ) return -1;
OrderedParticles testmode(outin.begin(), outin.end());
OrderedParticles::const_iterator dit = testmode.begin();
string testtag(in->name() + "->");
for( unsigned int i = 1; dit != testmode.end(); ++dit, ++i) {
testtag += (**dit).name();
if( i != 3 ) testtag += string(",");
}
if( testtag == refTag_ ) {
cc = false;
return 0;
}
else if ( testtag == refTagCC_ ) {
cc = true;
return 0;
}
else return -1;
}
bool GeneralThreeBodyDecayer::setDecayInfo(PDPtr incoming,
vector<PDPtr> outgoing,
const vector<TBDiagram> & process,
double symfac) {
// set the member variables from the info supplied
incoming_ = incoming;
- outgoing_ = outgoing;
+ outgoing_ = {outgoing[0],outgoing[1],outgoing[2]};
diagrams_ = process;
assert( outgoing_.size() == 3 );
// Construct reference tags for testing in modeNumber function
OrderedParticles refmode(outgoing_.begin(), outgoing_.end());
OrderedParticles::const_iterator dit = refmode.begin();
refTag_ = incoming_->name() + "->";
for( unsigned int i = 1; dit != refmode.end(); ++dit, ++i) {
refTag_ += (**dit).name();
if( i != 3 ) refTag_ += string(",");
}
//CC-mode
refmode.clear();
refTagCC_ = incoming_->CC() ? incoming_->CC()->name() :
incoming_->name();
refTagCC_ += "->";
for( unsigned int i = 0; i < 3; ++i ) {
if( outgoing_[i]->CC() ) refmode.insert( outgoing_[i]->CC() );
else refmode.insert( outgoing_[i] );
}
dit = refmode.begin();
for( unsigned int i = 1; dit != refmode.end(); ++dit , ++i) {
refTagCC_ += (**dit).name();
if( i != 3 ) refTagCC_ += string(",");
}
// check if intermeidates or only four point diagrams
bool intermediates(false);
for(auto diagram : diagrams_) {
if(diagram.intermediate) {
intermediates=true;
break;
}
}
if(!intermediates) {
incoming_= PDPtr();
outgoing_.clear();
generator()->log() << "Only four body diagrams for decay "
<< refTag_ << " in GeneralThreeBodyDecayer::"
<< "setDecayInfo(), omitting decay\n";
return false;
}
// set the colour factors and return the answer
if(setColourFactors(symfac)) return true;
incoming_= PDPtr();
outgoing_.clear();
return false;
}
void GeneralThreeBodyDecayer::doinit() {
- DecayIntegrator::doinit();
+ DecayIntegrator2::doinit();
if(outgoing_.empty()) return;
// create the phase space integrator
- tPDVector extpart(1,incoming_);
- extpart.insert(extpart.end(),outgoing_.begin(),outgoing_.end());
// create the integration channels for the decay
- DecayPhaseSpaceModePtr mode(new_ptr(DecayPhaseSpaceMode(extpart,this,true)));
- DecayPhaseSpaceChannelPtr newchannel;
+ PhaseSpaceModePtr mode(new_ptr(PhaseSpaceMode(incoming_,outgoing_,1.)));
// create the phase-space channels for the integration
unsigned int nmode(0);
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
if(diagrams_[ix].channelType==TBDiagram::fourPoint||
diagrams_[ix].channelType==TBDiagram::UNDEFINED) continue;
// create the new channel
- newchannel=new_ptr(DecayPhaseSpaceChannel(mode));
+ PhaseSpaceChannel newChannel(mode);
int jac = 0;
double power = 0.0;
if ( diagrams_[ix].intermediate->mass() == ZERO ||
diagrams_[ix].intermediate->width() == ZERO ) {
jac = 1;
power = -2.0;
}
- if(diagrams_[ix].channelType==TBDiagram::channel23) {
- newchannel->addIntermediate(extpart[0],0,0.0,-1,1);
- newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 2,3);
- }
- else if(diagrams_[ix].channelType==TBDiagram::channel13) {
- newchannel->addIntermediate(extpart[0],0,0.0,-1,2);
- newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 1,3);
- }
- else if(diagrams_[ix].channelType==TBDiagram::channel12) {
- newchannel->addIntermediate(extpart[0],0,0.0,-1,3);
- newchannel->addIntermediate(diagrams_[ix].intermediate,jac,power, 1,2);
- }
+ if(diagrams_[ix].channelType==TBDiagram::channel23)
+ newChannel = (newChannel,0,diagrams_[ix].intermediate,0,1,1,2,1,3);
+ else if(diagrams_[ix].channelType==TBDiagram::channel13)
+ newChannel = (newChannel,0,diagrams_[ix].intermediate,0,2,1,1,1,3);
+ else if(diagrams_[ix].channelType==TBDiagram::channel12)
+ newChannel = (newChannel,0,diagrams_[ix].intermediate,0,3,1,1,1,2);
+ if(jac==1)
+ newChannel.setJacobian(1,PhaseSpaceChannel::PhaseSpaceResonance::Power,power);
diagmap_.push_back(ix);
- mode->addChannel(newchannel);
+ mode->addChannel(newChannel);
++nmode;
}
if(nmode==0) {
- string mode = extpart[0]->PDGName() + "->";
- for(unsigned int ix=1;ix<extpart.size();++ix) mode += extpart[ix]->PDGName() + " ";
+ string mode = incoming_->PDGName() + "->";
+ for(unsigned int ix=0;ix<outgoing_.size();++ix) mode += outgoing_[ix]->PDGName() + " ";
throw Exception() << "No decay channels in GeneralThreeBodyDecayer::"
<< "doinit() for " << mode << "\n" << Exception::runerror;
}
// add the mode
- vector<double> wgt(nmode,1./double(nmode));
- addMode(mode,1.,wgt);
+ addMode(mode);
}
double GeneralThreeBodyDecayer::
threeBodyMatrixElement(const int imode, const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy m1,
const Energy m2, const Energy m3) const {
// calculate the momenta of the outgoing particles
Energy m0=sqrt(q2);
// energies
Energy eout[3] = {0.5*(q2+sqr(m1)-s1)/m0,
0.5*(q2+sqr(m2)-s2)/m0,
0.5*(q2+sqr(m3)-s3)/m0};
// magnitudes of the momenta
Energy pout[3] = {sqrt(sqr(eout[0])-sqr(m1)),
sqrt(sqr(eout[1])-sqr(m2)),
sqrt(sqr(eout[2])-sqr(m3))};
double cos2 = 0.5*(sqr(pout[0])+sqr(pout[1])-sqr(pout[2]))/pout[0]/pout[1];
double cos3 = 0.5*(sqr(pout[0])-sqr(pout[1])+sqr(pout[2]))/pout[0]/pout[2];
double sin2 = sqrt(1.-sqr(cos2)), sin3 = sqrt(1.-sqr(cos3));
- Lorentz5Momentum out[3]=
+ vector<Lorentz5Momentum> out =
{Lorentz5Momentum( ZERO , ZERO , pout[0] , eout[0] , m1),
Lorentz5Momentum( pout[1]*sin2 , ZERO , -pout[1]*cos2 , eout[1] , m2),
Lorentz5Momentum( -pout[2]*sin3 , ZERO , -pout[2]*cos3 , eout[2] , m3)};
- // create the incoming
- PPtr inpart=mode(imode)->externalParticles(0)->
+ // create the incoming particle
+ PPtr inpart=mode(imode)->incoming().first->
produceParticle(Lorentz5Momentum(sqrt(q2)));
- // and outgoing particles
- ParticleVector decay;
- for(unsigned int ix=1;ix<4;++ix)
- decay.push_back(mode(imode)->externalParticles(ix)->produceParticle(out[ix-1]));
// return the matrix element
- return me2(-1,*inpart,decay,Initialize);
+ return me2(-1,*inpart,outgoing_,out,Initialize);
}
double GeneralThreeBodyDecayer::brat(const DecayMode &, const Particle & p,
double oldbrat) const {
ParticleVector children = p.children();
if( children.size() != 3 || !p.data().widthGenerator() )
return oldbrat;
// partial width for this mode
Energy scale = p.mass();
Energy pwidth =
partialWidth( make_pair(p.dataPtr(), scale),
make_pair(children[0]->dataPtr(), children[0]->mass()),
make_pair(children[1]->dataPtr(), children[1]->mass()),
make_pair(children[2]->dataPtr(), children[2]->mass()) );
Energy width = p.data().widthGenerator()->width(p.data(), scale);
return pwidth/width;
}
Energy GeneralThreeBodyDecayer::partialWidth(PMPair inpart, PMPair outa,
PMPair outb, PMPair outc) const {
if(inpart.second<outa.second+outb.second+outc.second) return ZERO;
// create the object to calculate the width if it doesn't all ready exist
if(!widthCalc_) {
string tag = incoming_->name() + "->";
tag += outgoing_[0]->name() + "," + outgoing_[1]->name() + ","
+ outgoing_[2]->name() + ";";
DMPtr dm = generator()->findDecayMode(tag);
widthCalc_ = threeBodyMEIntegrator(*dm);
}
return widthCalc_->partialWidth(sqr(inpart.second));
}
void GeneralThreeBodyDecayer::
colourConnections(const Particle & parent,
const ParticleVector & out) const {
// first extract the outgoing particles and intermediate
PPtr inter;
ParticleVector outgoing;
if(!generateIntermediates()) {
outgoing=out;
}
else {
// find the diagram
unsigned int idiag = diagramMap()[mode(imode())->selectedChannel()];
PPtr child;
for(unsigned int ix=0;ix<out.size();++ix) {
if(out[ix]->children().empty()) child = out[ix];
else inter = out[ix];
}
outgoing.resize(3);
switch(diagrams_[idiag].channelType) {
case TBDiagram::channel23:
outgoing[0] = child;
outgoing[1] = inter->children()[0];
outgoing[2] = inter->children()[1];
break;
case TBDiagram::channel13:
outgoing[0] = inter->children()[0];
outgoing[1] = child;
outgoing[2] = inter->children()[1];
break;
case TBDiagram::channel12:
outgoing[0] = inter->children()[0];
outgoing[1] = inter->children()[1];
outgoing[2] = child;
break;
default:
throw Exception() << "unknown diagram type in GeneralThreeBodyDecayer::"
<< "colourConnections()" << Exception::runerror;
}
}
// extract colour of the incoming and outgoing particles
PDT::Colour inColour(parent.data().iColour());
vector<PDT::Colour> outColour;
vector<int> singlet,octet,triplet,antitriplet;
for(unsigned int ix=0;ix<outgoing.size();++ix) {
outColour.push_back(outgoing[ix]->data().iColour());
switch(outColour.back()) {
case PDT::Colour0 :
singlet.push_back(ix);
break;
case PDT::Colour3 :
triplet.push_back(ix);
break;
case PDT::Colour3bar:
antitriplet.push_back(ix);
break;
case PDT::Colour8 :
octet.push_back(ix);
break;
default:
throw Exception() << "Unknown colour for particle in GeneralThreeBodyDecayer::"
<< "colourConnections()" << Exception::runerror;
}
}
// colour neutral decaying particle
if ( inColour == PDT::Colour0) {
// options are all neutral or triplet/antitriplet+ neutral
if(singlet.size()==3) return;
else if(singlet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
// add intermediate if needed
if(inter&&inter->coloured()) {
if(inter->dataPtr()->iColour()==PDT::Colour3)
outgoing[triplet[0]]->colourLine()->addColoured(inter);
else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
}
}
else if(octet.size()==1&&triplet.size()==1&&antitriplet.size()==1) {
outgoing[ triplet[0]]->antiColourNeighbour(outgoing[octet[0]]);
outgoing[antitriplet[0]]-> colourNeighbour(outgoing[octet[0]]);
if(inter&&inter->coloured()) {
if(inter->dataPtr()->iColour()==PDT::Colour3)
outgoing[antitriplet[0]]->antiColourLine()->addColoured(inter);
else if(inter->dataPtr()->iColour()==PDT::Colour3bar)
outgoing[ triplet[0]]-> colourLine()->addAntiColoured(inter);
else if(inter->dataPtr()->iColour()==PDT::Colour8) {
outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
outgoing[ triplet[0]]-> colourLine()->addColoured(inter);
}
}
}
else if(triplet.size()==3) {
tColinePtr col[3] = {ColourLine::create(outgoing[0]),
ColourLine::create(outgoing[1]),
ColourLine::create(outgoing[2])};
col[0]->setSourceNeighbours(col[1],col[2]);
}
else if(antitriplet.size()==3) {
tColinePtr col[3] = {ColourLine::create(outgoing[0],true),
ColourLine::create(outgoing[1],true),
ColourLine::create(outgoing[2],true)};
col[0]->setSinkNeighbours(col[1],col[2]);
}
else {
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception()
<< "Unknown colour structure in GeneralThreeBodyDecayer::"
<< "colourConnections() for singlet decaying particle "
<< mode << Exception::runerror;
}
}
// colour triplet decaying particle
else if( inColour == PDT::Colour3) {
if(singlet.size()==2&&triplet.size()==1) {
outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
if(inter&&inter->coloured())
outgoing[triplet[0]]->colourLine()->addColoured(inter);
}
else if(antitriplet.size()==1&&triplet.size()==2) {
if(colourFlow()==0) {
outgoing[triplet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[1]]);
if(inter&&inter->coloured()) {
switch (inter->dataPtr()->iColour()) {
case PDT::Colour8:
inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[triplet[1]]->colourLine()->addAntiColoured(inter);
break;
default:
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception() << "Unknown colour for intermediate in "
<< "GeneralThreeBodyDecayer::"
<< "colourConnections() for "
<< "decaying colour triplet "
<< mode << Exception::runerror;
}
}
}
else {
outgoing[triplet[1]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[antitriplet[0]]->colourNeighbour(outgoing[triplet[0]]);
if(inter&&inter->coloured()) {
switch (inter->dataPtr()->iColour()) {
case PDT::Colour8:
inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[triplet[0]]->colourLine()->addAntiColoured(inter);
break;
default:
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception() << "Unknown colour for intermediate in "
<< "GeneralThreeBodyDecayer::"
<< "colourConnections() for "
<< "decaying colour triplet "
<< mode << Exception::runerror;
}
}
}
}
else if (singlet.size()==1&&triplet.size()==1&&octet.size()==1) {
if(inter) {
if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 ||
inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[octet[0]]->incomingColour(inter);
outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
}
else {
outgoing[octet[0]]->incomingColour(inter);
outgoing[octet[0]]->colourNeighbour(inter);
outgoing[triplet[0]]->incomingColour(inter);
}
}
else {
outgoing[octet[0]]->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[octet[0]]->colourNeighbour(outgoing[triplet[0]]);
}
}
else if (singlet.size()==1&&antitriplet.size()==2) {
tColinePtr col[2] = {ColourLine::create(outgoing[antitriplet[0]],true),
ColourLine::create(outgoing[antitriplet[1]],true)};
parent.colourLine()->setSinkNeighbours(col[0],col[1]);
}
else {
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception()
<< "Unknown colour structure in GeneralThreeBodyDecayer::"
<< "colourConnections() for triplet decaying particle "
<< mode << Exception::runerror;
}
}
else if( inColour == PDT::Colour3bar) {
if(singlet.size()==2&&antitriplet.size()==1) {
outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
else if(antitriplet.size()==2&&triplet.size()==1) {
if(colourFlow()==0) {
outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[1]]);
if(inter&&inter->coloured()) {
switch (inter->dataPtr()->iColour()) {
case PDT::Colour8:
inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
outgoing[antitriplet[1]]->antiColourLine()->addAntiColoured(inter);
break;
default:
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception() << "Unknown colour for intermediate in"
<< " GeneralThreeBodyDecayer::"
<< "colourConnections() for "
<< "decaying colour antitriplet "
<< mode << Exception::runerror;
}
}
}
else {
outgoing[antitriplet[1]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
outgoing[triplet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
if(inter&&inter->coloured()) {
switch (inter->dataPtr()->iColour()) {
case PDT::Colour8:
inter->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
break;
default:
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception() << "Unknown colour for intermediate in "
<< "GeneralThreeBodyDecayer::"
<< "colourConnections() for "
<< "decaying colour antitriplet "
<< mode << Exception::runerror;
}
}
}
}
else if (singlet.size()==1&&antitriplet.size()==1&&octet.size()==1) {
if(inter) {
if(inter->children()[0]->dataPtr()->iColour()==PDT::Colour8 ||
inter->children()[1]->dataPtr()->iColour()==PDT::Colour8) {
inter->incomingColour(const_ptr_cast<tPPtr>(&parent));
outgoing[octet[0]]->incomingAntiColour(inter);
outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
}
else {
outgoing[octet[0]]->incomingAntiColour(inter);
outgoing[octet[0]]->antiColourNeighbour(inter);
outgoing[antitriplet[0]]->incomingAntiColour(inter);
}
}
else {
outgoing[octet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
outgoing[octet[0]]->antiColourNeighbour(outgoing[antitriplet[0]]);
}
}
else if (singlet.size()==1&&triplet.size()==2) {
tColinePtr col[2] = {ColourLine::create(outgoing[triplet[0]]),
ColourLine::create(outgoing[triplet[1]])};
parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
}
else {
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception()
<< "Unknown colour structure in GeneralThreeBodyDecayer::"
<< "colourConnections() for anti-triplet decaying particle"
<< mode << Exception::runerror;
}
}
else if( inColour == PDT::Colour8) {
if(triplet.size()==1&&antitriplet.size()==1&&singlet.size()==1) {
outgoing[ triplet[0]]->incomingColour (const_ptr_cast<tPPtr>(&parent));
outgoing[antitriplet[0]]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
if(inter&&inter->coloured()) {
switch (inter->dataPtr()->iColour()) {
case PDT::Colour3:
outgoing[triplet[0]]->colourLine()->addColoured(inter);
break;
case PDT::Colour3bar:
outgoing[antitriplet[0]]->antiColourLine()->addAntiColoured(inter);
break;
default:
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception() << "Unknown colour for intermediate"
<< " in GeneralThreeBodyDecayer::"
<< "colourConnections() for "
<< "decaying colour octet "
<< mode << Exception::runerror;
}
}
}
else if(triplet.size()==3) {
tColinePtr col[2];
if(colourFlow()==0) {
outgoing[0]->incomingColour (const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[1]);
col[1] = ColourLine::create(outgoing[2]);
}
else if(colourFlow()==1) {
outgoing[1]->incomingColour (const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[0]);
col[1] = ColourLine::create(outgoing[2]);
}
else if(colourFlow()==2) {
outgoing[2]->incomingColour (const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[0]);
col[1] = ColourLine::create(outgoing[1]);
}
else
assert(false);
parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
}
else if(antitriplet.size()==3) {
tColinePtr col[2];
if(colourFlow()==0) {
outgoing[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[1],true);
col[1] = ColourLine::create(outgoing[2],true);
}
else if(colourFlow()==1) {
outgoing[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[0],true);
col[1] = ColourLine::create(outgoing[2],true);
}
else if(colourFlow()==2) {
outgoing[2]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
col[0] = ColourLine::create(outgoing[0],true);
col[1] = ColourLine::create(outgoing[1],true);
}
else
assert(false);
parent.colourLine()->setSinkNeighbours(col[0],col[1]);
}
else {
string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+ out[1]->PDGName() + " " + out[2]->PDGName();
throw Exception()
<< "Unknown colour structure in GeneralThreeBodyDecayer::"
<< "colourConnections() for octet decaying particle"
<< mode << Exception::runerror;
}
}
}
void GeneralThreeBodyDecayer::
constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
vector<Energy> & inwidth, vector<double> & inpow,
vector<double> & inweights) const {
// check if any intermediate photons
bool hasPhoton=false;
for(unsigned int iy=0;iy<diagmap_.size();++iy) {
unsigned int ix=diagmap_[iy];
if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma)
hasPhoton = true;
}
// loop over channels
Energy min = incoming()->mass();
int nchannel(0);
pair<int,Energy> imin[4]={make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV),
make_pair(-1,-1.*GeV),make_pair(-1,-1.*GeV)};
Energy absmin = -1e20*GeV;
int minType = -1;
for(unsigned int iy=0;iy<diagmap_.size();++iy) {
unsigned int ix=diagmap_[iy];
if(getProcessInfo()[ix].channelType==TBDiagram::fourPoint) continue;
Energy dm1(min-getProcessInfo()[ix].intermediate->mass());
Energy dm2(getProcessInfo()[ix].intermediate->mass());
int itype(0);
if (getProcessInfo()[ix].channelType==TBDiagram::channel23) {
dm1 -= outgoing()[0]->mass();
dm2 -= outgoing()[1]->mass()+outgoing()[2]->mass();
itype = 3;
}
else if(getProcessInfo()[ix].channelType==TBDiagram::channel13) {
dm1 -= outgoing()[1]->mass();
dm2 -= outgoing()[0]->mass()+outgoing()[2]->mass();
itype = 2;
}
else if(getProcessInfo()[ix].channelType==TBDiagram::channel12) {
dm1 -= outgoing()[2]->mass();
dm2 -= outgoing()[0]->mass()+outgoing()[1]->mass();
itype = 1;
}
if((dm1<ZERO||dm2<ZERO)&&!hasPhoton) {
if (imin[itype].first < 0 ||
(dm1<ZERO && imin[itype].second < dm1) ) {
imin[itype] = make_pair(ix,dm1);
if(dm1<ZERO&&absmin<dm1) {
absmin = dm1;
minType = itype;
}
}
continue;
}
if(getProcessInfo()[ix].intermediate->id()!=ParticleID::gamma) {
intype.push_back(itype);
inpow.push_back(0.);
inmass.push_back(getProcessInfo()[ix].intermediate->mass());
inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[ix].intermediate->width());
++nchannel;
}
else if(getProcessInfo()[ix].intermediate->id()==ParticleID::gamma) {
intype.push_back(itype);
inpow.push_back(-2.);
inmass.push_back(-1.*GeV);
inwidth.push_back(-1.*GeV);
++nchannel;
}
}
// physical poles, use them and return
if(nchannel>0) {
inweights = vector<double>(nchannel,1./double(nchannel));
return;
}
// use shallowest pole
else if(intOpt_==1&&minType>0&&getProcessInfo()[imin[minType].first].intermediate->id()!=ParticleID::gamma) {
intype.push_back(minType);
inpow.push_back(0.);
inmass.push_back(getProcessInfo()[imin[minType].first].intermediate->mass());
inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[minType].first].intermediate->width());
inweights = vector<double>(1,1.);
return;
}
for(unsigned int ix=1;ix<4;++ix) {
if(imin[ix].first>=0) {
intype.push_back(ix);
if(getProcessInfo()[imin[ix].first].intermediate->id()!=ParticleID::gamma) {
inpow.push_back(0.);
inmass.push_back(getProcessInfo()[imin[ix].first].intermediate->mass());
inwidth.push_back(widthOption() ==3 ? ZERO : getProcessInfo()[imin[ix].first].intermediate->width());
}
else {
inpow.push_back(-2.);
inmass.push_back(-1.*GeV);
inwidth.push_back(-1.*GeV);
}
++nchannel;
}
}
inweights = vector<double>(nchannel,1./double(nchannel));
}
bool GeneralThreeBodyDecayer::setColourFactors(double symfac) {
string name = incoming_->PDGName() + "->";
vector<int> sng,trip,atrip,oct;
unsigned int iloc(0);
- for(vector<PDPtr>::const_iterator it = outgoing_.begin();
+ for(tPDVector::const_iterator it = outgoing_.begin();
it != outgoing_.end();++it) {
name += (**it).PDGName() + " ";
if ((**it).iColour() == PDT::Colour0 ) sng.push_back(iloc) ;
else if((**it).iColour() == PDT::Colour3 ) trip.push_back(iloc) ;
else if((**it).iColour() == PDT::Colour3bar ) atrip.push_back(iloc);
else if((**it).iColour() == PDT::Colour8 ) oct.push_back(iloc) ;
++iloc;
}
// colour neutral decaying particle
if ( incoming_->iColour() == PDT::Colour0) {
// options are all neutral or triplet/antitriplet+ neutral
if(sng.size()==3) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,1.));
colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
}
else if(sng.size()==1&&trip.size()==1&&atrip.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,3.));
colourLargeNC_ = vector<DVector>(1,DVector(1,3.));
}
else if(trip.size()==1&&atrip.size()==1&&oct.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,4.));
colourLargeNC_ = vector<DVector>(1,DVector(1,4.));
}
else if( trip.size() == 3 || atrip.size() == 3 ) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,6.));
colourLargeNC_ = vector<DVector>(1,DVector(1,6.));
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
double sign = diagrams_[ix].channelType == TBDiagram::channel13 ? -1. : 1.;
diagrams_[ix]. colourFlow = vector<CFPair>(1,make_pair(1,sign));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,sign));
}
}
else {
generator()->log() << "Unknown colour flow structure for "
<< "colour neutral decay "
<< name << " in GeneralThreeBodyDecayer::"
<< "setColourFactors(), omitting decay\n";
return false;
}
}
// colour triplet decaying particle
else if( incoming_->iColour() == PDT::Colour3) {
if(sng.size()==2&&trip.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,1.));
colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
}
else if(trip.size()==2&&atrip.size()==1) {
nflow_ = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 3.; colour_[0][1] = 1.;
colour_[1][0] = 1.; colour_[1][1] = 3.;
colourLargeNC_.clear();
colourLargeNC_.resize(2,DVector(2,0.));
colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.;
// sort out the contribution of the different diagrams to the colour
// flows
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
// colour singlet intermediate
if(diagrams_[ix].intermediate->iColour()==PDT::Colour0) {
if(diagrams_[ix].channelType==trip[0]) {
diagrams_[ix]. colourFlow = vector<CFPair>(1,make_pair(1,1.));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
}
else {
diagrams_[ix].colourFlow = vector<CFPair>(1,make_pair(2,1.));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
}
}
// colour octet intermediate
else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) {
if(diagrams_[ix].channelType==trip[0]) {
vector<CFPair> flow(1,make_pair(2, 0.5 ));
diagrams_[ix].largeNcColourFlow = flow;
flow.push_back( make_pair(1,-1./6.));
diagrams_[ix].colourFlow=flow;
}
else {
vector<CFPair> flow(1,make_pair(1, 0.5 ));
diagrams_[ix].largeNcColourFlow = flow;
flow.push_back( make_pair(2,-1./6.));
diagrams_[ix].colourFlow=flow;
}
}
else {
generator()->log() << "Unknown colour for the intermediate in "
<< "triplet -> triplet triplet antitriplet in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
}
}
else if(trip.size()==1&&oct.size()==1&&sng.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,4./3.));
colourLargeNC_ = vector<DVector>(1,DVector(1,4./3.));
}
else if(sng.size()==1&&atrip.size()==2) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,2.));
colourLargeNC_ = vector<DVector>(1,DVector(1,2.));
}
else {
generator()->log() << "Unknown colour structure for "
<< "triplet decay in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
}
// colour antitriplet decaying particle
else if( incoming_->iColour() == PDT::Colour3bar) {
if(sng.size()==2&&atrip.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,1.));
colourLargeNC_ = vector<DVector>(1,DVector(1,1.));
}
else if(atrip.size()==2&&trip.size()==1) {
nflow_ = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = 3.; colour_[0][1] = 1.;
colour_[1][0] = 1.; colour_[1][1] = 3.;
colourLargeNC_.clear();
colourLargeNC_.resize(2,DVector(2,0.));
colourLargeNC_[0][0] = 3.; colourLargeNC_[1][1] = 3.;
// sort out the contribution of the different diagrams to the colour
// flows
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
// colour singlet intermediate
if(diagrams_[ix].intermediate->iColour()==PDT::Colour0) {
if(diagrams_[ix].channelType==atrip[0]) {
diagrams_[ix]. colourFlow = vector<CFPair>(1,make_pair(1,1.));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
}
else {
diagrams_[ix].colourFlow = vector<CFPair>(1,make_pair(2,1.));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
}
}
// colour octet intermediate
else if(diagrams_[ix].intermediate->iColour()==PDT::Colour8) {
if(diagrams_[ix].channelType==atrip[0]) {
vector<CFPair> flow(1,make_pair(2, 0.5 ));
diagrams_[ix].largeNcColourFlow = flow;
flow.push_back( make_pair(1,-1./6.));
diagrams_[ix].colourFlow=flow;
}
else {
vector<CFPair> flow(1,make_pair(1, 0.5 ));
diagrams_[ix].largeNcColourFlow = flow;
flow.push_back( make_pair(2,-1./6.));
diagrams_[ix].colourFlow=flow;
}
}
else {
generator()->log() << "Unknown colour for the intermediate in "
<< "antitriplet -> antitriplet antitriplet triplet in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
}
}
else if(atrip.size()==1&&oct.size()==1&&sng.size()==1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,4./3.));
colourLargeNC_ = vector<DVector>(1,DVector(1,4./3.));
}
else if(sng.size()==1&&trip.size()==2) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,2.));
colourLargeNC_ = vector<DVector>(1,DVector(1,2.));
}
else {
generator()->log() << "Unknown colour antitriplet decay in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
}
// colour octet particle
else if( incoming_->iColour() == PDT::Colour8) {
// triplet antitriplet
if(trip.size() == 1 && atrip.size() == 1 && sng.size() == 1) {
nflow_ = 1;
colour_ = vector<DVector>(1,DVector(1,0.5));
colourLargeNC_ = vector<DVector>(1,DVector(1,0.5));
}
// three (anti)triplets
else if(trip.size()==3||atrip.size()==3) {
nflow_ = 3;
colour_ = vector<DVector>(3,DVector(3,0.));
colourLargeNC_ = vector<DVector>(3,DVector(3,0.));
colour_[0][0] = 1.; colour_[1][1] = 1.; colour_[2][2] = 1.;
colour_[0][1] = -0.5; colour_[1][0] = -0.5;
colour_[0][2] = -0.5; colour_[2][0] = -0.5;
colour_[1][2] = -0.5; colour_[2][1] = -0.5;
colourLargeNC_ = vector<DVector>(3,DVector(3,0.));
colourLargeNC_[0][0] = 1.; colourLargeNC_[1][1] = 1.; colourLargeNC_[2][2] = 1.;
// sett the factors for the diagrams
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
tPDPtr inter = diagrams_[ix].intermediate;
if(inter->CC()) inter = inter->CC();
unsigned int io[2]={1,2};
double sign = diagrams_[ix].channelType == TBDiagram::channel13 ? -1. : 1.;
for(unsigned int iy=0;iy<3;++iy) {
if (iy==1) io[0]=0;
else if(iy==2) io[1]=1;
tPDVector decaylist = diagrams_[ix].vertices.second->search(iy, inter);
if(decaylist.empty()) continue;
bool found=false;
for(unsigned int iz=0;iz<decaylist.size();iz+=3) {
if(decaylist[iz+io[0]]->id()==diagrams_[ix].outgoingPair.first &&
decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.second) {
sign *= 1.;
found = true;
}
else if(decaylist[iz+io[0]]->id()==diagrams_[ix].outgoingPair.second &&
decaylist[iz+io[1]]->id()==diagrams_[ix].outgoingPair.first ) {
sign *= -1.;
found = true;
}
}
if(found) {
if(iy==1) sign *=-1.;
break;
}
}
diagrams_[ix]. colourFlow = vector<CFPair>(1,make_pair(diagrams_[ix].channelType+1,sign));
diagrams_[ix].largeNcColourFlow = vector<CFPair>(1,make_pair(diagrams_[ix].channelType+1,sign));
}
}
// unknown
else {
generator()->log() << "Unknown colour octet decay in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
}
else if (incoming_->iColour() == PDT::Colour6 ) {
generator()->log() << "Unknown colour sextet decay in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
else if (incoming_->iColour() == PDT::Colour6bar ) {
generator()->log() << "Unknown colour anti-sextet decay in "
<< "GeneralThreeBodyDecayer::setColourFactors()"
<< " for " << name << " omitting decay\n";
return false;
}
assert(nflow_ != 999);
for(unsigned int ix=0;ix<nflow_;++ix) {
for(unsigned int iy=0;iy<nflow_;++iy) {
colour_ [ix][iy] /= symfac;
colourLargeNC_[ix][iy] /= symfac;
}
}
if( Debug::level > 1 ) {
generator()->log() << "Mode: " << name << " has colour factors\n";
for(unsigned int ix=0;ix<nflow_;++ix) {
for(unsigned int iy=0;iy<nflow_;++iy) {
generator()->log() << colour_[ix][iy] << " ";
}
generator()->log() << "\n";
}
for(unsigned int ix=0;ix<diagrams_.size();++ix) {
generator()->log() << "colour flow for diagram : " << ix;
for(unsigned int iy=0;iy<diagrams_[ix].colourFlow.size();++iy)
generator()->log() << "(" << diagrams_[ix].colourFlow[iy].first << ","
<< diagrams_[ix].colourFlow[iy].second << "); ";
generator()->log() << "\n";
}
}
return true;
}
diff --git a/Decay/General/GeneralThreeBodyDecayer.h b/Decay/General/GeneralThreeBodyDecayer.h
--- a/Decay/General/GeneralThreeBodyDecayer.h
+++ b/Decay/General/GeneralThreeBodyDecayer.h
@@ -1,321 +1,321 @@
// -*- C++ -*-
#ifndef HERWIG_GeneralThreeBodyDecayer_H
#define HERWIG_GeneralThreeBodyDecayer_H
//
// This is the declaration of the GeneralThreeBodyDecayer class.
//
-#include "Herwig/Decay/DecayIntegrator.h"
+#include "Herwig/Decay/DecayIntegrator2.h"
#include "Herwig/Models/General/TBDiagram.h"
#include "GeneralThreeBodyDecayer.fh"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the GeneralThreeBodyDecayer class.
*
* @see \ref GeneralThreeBodyDecayerInterfaces "The interfaces"
* defined for GeneralThreeBodyDecayer.
*/
-class GeneralThreeBodyDecayer: public DecayIntegrator {
+class GeneralThreeBodyDecayer: public DecayIntegrator2 {
public:
/** A ParticleData ptr and (possible) mass pair.*/
typedef pair<tcPDPtr, Energy> PMPair;
public:
/**
* The default constructor.
*/
GeneralThreeBodyDecayer() : nflow_(999), widthOpt_(1),
refTag_(), refTagCC_(), iflow_(999),
intOpt_(0), relerr_(1e-2)
{}
/** @name Virtual functions required by the Decayer class. */
//@{
/**
* For a given decay mode and a given particle instance, perform the
* decay and return the decay products. As this is the base class this
* is not implemented.
* @return The vector of particles produced in the decay.
*/
virtual ParticleVector decay(const Particle & parent,
const tPDVector & children) const;
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,const tPDVector & children) const;
/**
* The matrix element to be integrated for the three-body decays as a function
* of the invariant masses of pairs of the outgoing particles.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The matrix element
*/
virtual double threeBodyMatrixElement(const int imode, const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy m1,
const Energy m2, const Energy m3) const;
/**
* Function to return partial Width
* @param inpart The decaying particle.
* @param outa First decay product.
* @param outb Second decay product.
* @param outc Third decay product.
*/
virtual Energy partialWidth(PMPair inpart, PMPair outa,
PMPair outb, PMPair outc) const;
/**
* An overidden member to calculate a branching ratio for a certain
* particle instance.
* @param dm The DecayMode of the particle
* @param p The particle object
* @param oldbrat The branching fraction given in the DecayMode object
*/
virtual double brat(const DecayMode & dm, const Particle & p,
double oldbrat) const;
//@}
/**
* Set the diagrams
*/
bool setDecayInfo(PDPtr incoming,vector<PDPtr> outgoing,
const vector<TBDiagram> & process,
double symfac);
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 Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
protected:
/**
* Access the TBDiagrams that store the required information
* to create the diagrams
*/
const vector<TBDiagram> & getProcessInfo() const {
return diagrams_;
}
/**
* Incoming particle
*/
PDPtr incoming() const { return incoming_; }
/**
* Outgoing particles
*/
- const vector<PDPtr> & outgoing() const { return outgoing_; }
+ const vector<tPDPtr> & outgoing() const { return outgoing_; }
/**
* Number of colour flows
*/
unsigned int numberOfFlows() const { return nflow_; }
/**
* Set up the colour factors
*/
bool setColourFactors(double symfac);
/**
* Return the matrix of colour factors
*/
const vector<DVector> & getColourFactors() const { return colour_; }
/**
* Return the matrix of colour factors
*/
const vector<DVector> & getLargeNcColourFactors() const {
return colourLargeNC_;
}
/**
* Get the mapping between the phase-space channel and the diagram
*/
const vector<unsigned int> & diagramMap() const {
return diagmap_;
}
/**
* Option for the handling of the widths of the intermediate particles
*/
unsigned int widthOption() const { return widthOpt_; }
/**
* Set colour connections
* @param parent Parent particle
* @param out Particle vector containing particles to
* connect colour lines
*/
void colourConnections(const Particle & parent,
const ParticleVector & out) const;
/**
* Method to construct the channels for the integrator to give the partial width
* @param intype Types of the channels
* @param inmass Mass for the channels
* @param inwidth Width for the channels
* @param inpow Power for the channels
* @param inweights Weights for the channels
*/
void constructIntegratorChannels(vector<int> & intype, vector<Energy> & inmass,
vector<Energy> & inwidth, vector<double> & inpow,
vector<double> & inweights) const;
/**
* Set the colour flow
* @param flow The value for the colour flow
*/
void colourFlow(unsigned int flow) const { iflow_ = flow; }
/**
* Set the colour flow
*/
unsigned int const & colourFlow() const { return iflow_; }
/**
* Relative error for GQ integration
*/
double relativeError() const {return relerr_;}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
GeneralThreeBodyDecayer & operator=(const GeneralThreeBodyDecayer &);
private:
/**
* Store the incoming particle
*/
PDPtr incoming_;
/**
* Outgoing particles
*/
- vector<PDPtr> outgoing_;
+ vector<tPDPtr> outgoing_;
/**
* Store the diagrams for the decay
*/
vector<TBDiagram> diagrams_;
/**
* Map between the diagrams and the phase-space channels
*/
vector<unsigned int> diagmap_;
/**
* Store colour factors for ME calc.
*/
vector<DVector> colour_;
/**
* Store cololur factors for ME calc at large N_c
*/
vector<DVector> colourLargeNC_;
/**
* The number of colourflows.
*/
unsigned int nflow_;
/**
* Reference to object to calculate the partial width
*/
mutable WidthCalculatorBasePtr widthCalc_;
/**
* Option for the treatment of the widths
*/
unsigned int widthOpt_;
/**
* Store a decay tag for this mode that can be tested when
* trying to determine whether it can be generated by
* this Decayer
*/
string refTag_;
/**
* Store a decay tag for the cc-mode that can be tested when
* trying to determine whether it can be generated by
* this Decayer
*/
string refTagCC_;
/**
* The colour flow
*/
mutable unsigned int iflow_;
/**
* Option for the construction of the gaussian integrator
*/
unsigned int intOpt_;
/**
* Relative error for GQ integration of partial width
*/
double relerr_;
};
}
#endif /* HERWIG_GeneralThreeBodyDecayer_H */
diff --git a/Decay/General/SFFDecayer.h b/Decay/General/SFFDecayer.h
--- a/Decay/General/SFFDecayer.h
+++ b/Decay/General/SFFDecayer.h
@@ -1,241 +1,236 @@
// -*- 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 "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
* partialWidth().
*
* @see GeneralTwoBodyDecayer
*/
class SFFDecayer: public GeneralTwoBodyDecayer {
public:
-
- /**
- * The default constructor.
- */
- SFFDecayer() {}
/**
* Return the matrix element squared for a given mode and phase-space channel.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param outgoing The particles produced in the decay
* @param momenta The momenta of the particles produced in the decay
* @param meopt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
double me2(const int ichan,const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const;
/**
* Construct the SpinInfos for the particles produced in the decay
*/
virtual void constructSpinInfo(const Particle & part,
ParticleVector outgoing) const;
/**
* 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/StoFFVDecayer.cc b/Decay/General/StoFFVDecayer.cc
--- a/Decay/General/StoFFVDecayer.cc
+++ b/Decay/General/StoFFVDecayer.cc
@@ -1,389 +1,393 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StoFFVDecayer class.
//
#include "StoFFVDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG;
using namespace ThePEG::Helicity;
IBPtr StoFFVDecayer::clone() const {
return new_ptr(*this);
}
IBPtr StoFFVDecayer::fullclone() const {
return new_ptr(*this);
}
void StoFFVDecayer::persistentOutput(PersistentOStream & os) const {
os << sca_ << fer_ << vec_ << RSfer_ << four_;
}
void StoFFVDecayer::persistentInput(PersistentIStream & is, int) {
is >> sca_ >> fer_ >> vec_ >> RSfer_ >> four_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<StoFFVDecayer,GeneralThreeBodyDecayer>
describeHerwigStoFFVDecayer("Herwig::StoFFVDecayer", "Herwig.so");
void StoFFVDecayer::Init() {
static ClassDocumentation<StoFFVDecayer> documentation
("The StoFFVDecayer class implements the general decay of a scalar to "
"a two fermions and a vector.");
}
WidthCalculatorBasePtr StoFFVDecayer::
threeBodyMEIntegrator(const DecayMode & ) const {
vector<int> intype;
vector<Energy> inmass,inwidth;
vector<double> inpow,inweights;
constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights);
return new_ptr(ThreeBodyAllOnCalculator<StoFFVDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,0,
outgoing()[0]->mass(),outgoing()[1]->mass(),
outgoing()[2]->mass(),relativeError()));
}
void StoFFVDecayer::doinit() {
GeneralThreeBodyDecayer::doinit();
if(outgoing().empty()) return;
unsigned int ndiags = getProcessInfo().size();
sca_.resize(ndiags);
fer_.resize(ndiags);
RSfer_.resize(ndiags);
vec_.resize(ndiags);
four_.resize(ndiags);
for(unsigned int ix = 0;ix < ndiags; ++ix) {
TBDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
// four point vertex
if(!offshell) {
four_[ix] = dynamic_ptr_cast<AbstractFFVSVertexPtr>(current.vertices.first);
continue;
}
if( offshell->CC() ) offshell = offshell->CC();
if(offshell->iSpin() == PDT::Spin0) {
AbstractVSSVertexPtr vert1 = dynamic_ptr_cast<AbstractVSSVertexPtr>
(current.vertices.first);
AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in StoFFVDecayer::doinit()"
<< Exception::runerror;
sca_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1Half) {
AbstractFFSVertexPtr vert1 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a fermion diagram in StoFFVDecayer::doinit()"
<< Exception::runerror;
fer_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractVVSVertexPtr vert1 = dynamic_ptr_cast<AbstractVVSVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a vector diagram in StoFFVDecayer::doinit()"
<< Exception::runerror;
vec_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin3Half) {
AbstractRFSVertexPtr vert1 = dynamic_ptr_cast<AbstractRFSVertexPtr>
(current.vertices.first);
AbstractRFVVertexPtr vert2 = dynamic_ptr_cast<AbstractRFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a RS fermion diagram in StoFFVDecayer::doinit()"
<< Exception::runerror;
RSfer_[ix] = make_pair(vert1, vert2);
}
}
}
-double StoFFVDecayer::me2(const int ichan, const Particle & inpart,
- const ParticleVector & decay, MEOption meopt) const {
+void StoFFVDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ ScalarWaveFunction::
+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ for(unsigned int ix=0;ix<decay.size();++ix) {
+ if(decay[ix]->dataPtr()->iSpin()==PDT::Spin1) {
+ VectorWaveFunction::constructSpinInfo(outVector_,decay[ix],
+ Helicity::outgoing,true,false);
+ }
+ else {
+ SpinorWaveFunction::
+ constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
+ }
+ }
+}
+
+double StoFFVDecayer::me2(const int ichan, const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ MEOption meopt) const {
// particle or CC of particle
- bool cc = (*getProcessInfo().begin()).incoming != inpart.id();
+ bool cc = (*getProcessInfo().begin()).incoming != part.id();
// special handling or first/last call
if(meopt==Initialize) {
ScalarWaveFunction::
- calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming);
- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),
+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),
Helicity::incoming);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- ScalarWaveFunction::
- constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- for(unsigned int ix=0;ix<decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin()==PDT::Spin1) {
- VectorWaveFunction::constructSpinInfo(outVector_,decay[ix],
- Helicity::outgoing,true,false);
- }
- else {
- SpinorWaveFunction::
- constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
- }
- }
- }
unsigned int ivec(0);
bool massless(false);
- for(unsigned int ix = 0; ix < decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin() == PDT::Spin1) {
+ for(unsigned int ix = 0; ix < outgoing.size();++ix) {
+ if(outgoing[ix]->iSpin() == PDT::Spin1) {
ivec = ix;
- massless = decay[ivec]->mass()==ZERO;
+ massless = outgoing[ivec]->mass()==ZERO;
VectorWaveFunction::
- calculateWaveFunctions(outVector_, decay[ix], Helicity::outgoing,massless);
+ calculateWaveFunctions(outVector_, momenta[ix],outgoing[ix], Helicity::outgoing,massless);
}
else {
SpinorWaveFunction::
- calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing);
+ calculateWaveFunctions(outspin_[ix].first,momenta[ix],outgoing[ix],Helicity::outgoing);
outspin_[ix].second.resize(2);
// Need a ubar and a v spinor
if(outspin_[ix].first[0].wave().Type() == SpinorType::u) {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].first[iy].conjugate();
}
}
else {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].second[iy].conjugate();
}
}
}
}
const vector<vector<double> > cfactors(getColourFactors());
const vector<vector<double> > nfactors(getLargeNcColourFactors());
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
const size_t ncf(numberOfFlows());
vector<Complex> flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.));
// setup the DecayMatrixElement
vector<GeneralDecayMEPtr>
mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0,
ivec == 0 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 1 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 2 ? PDT::Spin1 : PDT::Spin1Half)));
vector<GeneralDecayMEPtr>
mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0,
ivec == 0 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 1 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 2 ? PDT::Spin1 : PDT::Spin1Half)));
//the channel possiblities
static const unsigned int out2[3] = {1,0,0}, out3[3] = {2,2,1};
for(unsigned int s1 = 0; s1 < 2; ++s1) {
for(unsigned int s2 = 0; s2 < 2; ++s2) {
for(unsigned int v1 = 0; v1 < 3; ++v1) {
if(massless&&v1==1) ++v1;
flows = vector<Complex>(ncf, Complex(0.));
largeflows = vector<Complex>(ncf, Complex(0.));
unsigned int idiag(0);
Complex diag;
for(vector<TBDiagram>::const_iterator dit=getProcessInfo().begin();
dit!=getProcessInfo().end();++dit) {
// channels if selecting
if( ichan >= 0 && diagramMap()[ichan] != idiag ) {
++idiag;
continue;
}
tcPDPtr offshell = dit->intermediate;
if(offshell) {
if(cc&&offshell->CC()) offshell=offshell->CC();
unsigned int o2(out2[dit->channelType]), o3(out3[dit->channelType]);
double sign = (o3 < o2) ? 1. : -1.;
// intermediate scalar
if(offshell->iSpin() == PDT::Spin0) {
ScalarWaveFunction inters = sca_[idiag].first->
evaluate(scale, widthOption(), offshell, outVector_[v1], swave_);
unsigned int h1(s1),h2(s2);
if(o2 > o3) swap(h1, h2);
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag = -sign*sca_[idiag].second->
evaluate(scale,outspin_[o2].first[h1],
outspin_[o3].second[h2],inters);
}
else {
diag = sign*sca_[idiag].second->
evaluate(scale, outspin_[o3].first [h2],
outspin_[o2].second[h1],inters);
}
}
// intermediate fermion
else if(offshell->iSpin() == PDT::Spin1Half) {
- int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half)
+ int iferm = (outgoing[o2]->iSpin() == PDT::Spin1Half)
? o2 : o3;
unsigned int h1(s1),h2(s2);
if(dit->channelType > iferm) swap(h1, h2);
sign = iferm < dit->channelType ? 1. : -1.;
- if((decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) ||
- (decay[dit->channelType]->id()*offshell->id()>0)) {
+ if((outgoing[dit->channelType]->id() < 0 && outgoing[iferm]->id() > 0 ) ||
+ (outgoing[dit->channelType]->id()*offshell->id()>0)) {
SpinorWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].first[h1], swave_);
diag = -sign*fer_[idiag].second->
evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]);
}
else {
SpinorBarWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].second[h1],swave_);
diag = sign*fer_[idiag].second->
evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]);
}
}
// intermediate vector
else if(offshell->iSpin() == PDT::Spin1) {
VectorWaveFunction interv = vec_[idiag].first->
evaluate(scale, widthOption(), offshell, outVector_[v1], swave_);
unsigned int h1(s1),h2(s2);
if(o2 > o3) swap(h1,h2);
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag =-sign*vec_[idiag].second->
evaluate(scale, outspin_[o2].first[h1],
outspin_[o3].second[h2], interv);
}
else {
diag = sign*vec_[idiag].second->
evaluate(scale, outspin_[o3].first[h2],
outspin_[o2].second[h1], interv);
}
}
// intermediate RS fermion
else if(offshell->iSpin() == PDT::Spin3Half) {
- int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half)
+ int iferm = (outgoing[o2]->iSpin() == PDT::Spin1Half)
? o2 : o3;
unsigned int h1(s1),h2(s2);
if(dit->channelType > iferm) swap(h1, h2);
sign = iferm < dit->channelType ? 1. : -1.;
- if((decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) ||
- (decay[dit->channelType]->id()*offshell->id()>0)) {
+ if((outgoing[dit->channelType]->id() < 0 && outgoing[iferm]->id() > 0 ) ||
+ (outgoing[dit->channelType]->id()*offshell->id()>0)) {
RSSpinorWaveFunction inters = RSfer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].first[h1], swave_);
diag = -sign*RSfer_[idiag].second->
evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]);
}
else {
RSSpinorBarWaveFunction inters = RSfer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].second[h1],swave_);
diag = sign*RSfer_[idiag].second->
evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]);
}
}
// unknown
else throw Exception()
<< "Unknown intermediate in StoFFVDecayer::me2()"
<< Exception::runerror;
}
else {
unsigned int o2 = ivec > 0 ? 0 : 1;
unsigned int o3 = ivec < 2 ? 2 : 1;
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag =-four_[idiag]->
evaluate(scale, outspin_[o2].first[s1],
outspin_[o3].second[s2], outVector_[v1], swave_);
}
else {
diag = four_[idiag]->
evaluate(scale, outspin_[o3].first[s2],
outspin_[o2].second[s1], outVector_[v1], swave_);
}
}
// matrix element for the different colour flows
if(ichan < 0) {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
else {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1 != colourFlow()) continue;
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
++idiag;
} //end of diagrams
// now add the flows to the me2 with appropriate colour factors
for(unsigned int ix = 0; ix < ncf; ++ix) {
if ( ivec == 0 ) {
(*mes[ix])(0, v1, s1, s2) = flows[ix];
(*mel[ix])(0, v1, s1, s2) = largeflows[ix];
}
else if( ivec == 1 ) {
(*mes[ix])(0, s1, v1, s2) = flows[ix];
(*mel[ix])(0, s1, v1, s2) = largeflows[ix];
}
else if( ivec == 2 ) {
(*mes[ix])(0, s1, s2, v1) = flows[ix];
(*mel[ix])(0, s1, s2, v1) = largeflows[ix];
}
}
}
}
}
double me2(0.);
if(ichan < 0) {
vector<double> pflows(ncf,0.);
for(unsigned int ix = 0; ix < ncf; ++ix) {
for(unsigned int iy = 0; iy < ncf; ++ iy) {
double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real();
me2 += con;
if(ix == iy) {
con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real();
pflows[ix] += con;
}
}
}
double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.));
ptotal *= UseRandom::rnd();
for(unsigned int ix = 0;ix < pflows.size(); ++ix) {
if(ptotal <= pflows[ix]) {
colourFlow(ix);
ME(mes[ix]);
break;
}
ptotal -= pflows[ix];
}
}
else {
unsigned int iflow = colourFlow();
me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real();
}
// return the matrix element squared
return me2;
}
diff --git a/Decay/General/StoFFVDecayer.h b/Decay/General/StoFFVDecayer.h
--- a/Decay/General/StoFFVDecayer.h
+++ b/Decay/General/StoFFVDecayer.h
@@ -1,162 +1,171 @@
// -*- C++ -*-
#ifndef THEPEG_StoFFVDecayer_H
#define THEPEG_StoFFVDecayer_H
//
// This is the declaration of the StoFFVDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVSVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the StoFFVDecayer class.
*
* @see \ref StoFFVDecayerInterfaces "The interfaces"
* defined for StoFFVDecayer.
*/
class StoFFVDecayer: public GeneralThreeBodyDecayer {
public:
-
+
/**
- * 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;
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of
* calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
StoFFVDecayer & operator=(const StoFFVDecayer &);
private:
/**
* Store the vertices for fermion intrermediate
*/
vector<pair<AbstractFFSVertexPtr, AbstractFFVVertexPtr> > fer_;
/**
* Store the vertices for fermion intrermediate
*/
vector<pair<AbstractRFSVertexPtr, AbstractRFVVertexPtr> > RSfer_;
/**
* Store the vertices for scalar intrermediate
*/
vector<pair<AbstractVSSVertexPtr, AbstractFFSVertexPtr> > sca_;
/**
* Store the vertices for vector intrermediate
*/
vector<pair<AbstractVVSVertexPtr, AbstractFFVVertexPtr> > vec_;
/**
* Store the vertices for 4-point diagrams
*/
vector<AbstractFFVSVertexPtr> four_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Scalar wavefunction
*/
mutable ScalarWaveFunction swave_;
/**
* Vector wavefunction
*/
mutable vector<VectorWaveFunction> outVector_;
/**
* Spinor wavefunctions
*/
mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > outspin_[3];
};
}
#endif /* THEPEG_StoFFVDecayer_H */
diff --git a/Decay/General/StoSFFDecayer.cc b/Decay/General/StoSFFDecayer.cc
--- a/Decay/General/StoSFFDecayer.cc
+++ b/Decay/General/StoSFFDecayer.cc
@@ -1,424 +1,425 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StoSFFDecayer class.
//
#include "StoSFFDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG;
using namespace ThePEG::Helicity;
IBPtr StoSFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr StoSFFDecayer::fullclone() const {
return new_ptr(*this);
}
void StoSFFDecayer::persistentOutput(PersistentOStream & os) const {
os << sca_ << fer_ << vec_ << ten_ << RSfer_ << four_;
}
void StoSFFDecayer::persistentInput(PersistentIStream & is, int) {
is >> sca_ >> fer_ >> vec_ >> ten_ >> RSfer_ >> four_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<StoSFFDecayer,GeneralThreeBodyDecayer>
describeHerwigStoSFFDecayer("Herwig::StoSFFDecayer", "Herwig.so");
void StoSFFDecayer::Init() {
static ClassDocumentation<StoSFFDecayer> documentation
("The StoSFFDecayer class implements the general decay of a scalar to "
"a scalar and two fermions.");
}
WidthCalculatorBasePtr StoSFFDecayer::
threeBodyMEIntegrator(const DecayMode & ) const {
vector<int> intype;
vector<Energy> inmass,inwidth;
vector<double> inpow,inweights;
constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights);
return new_ptr(ThreeBodyAllOnCalculator<StoSFFDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,0,
outgoing()[0]->mass(),outgoing()[1]->mass(),outgoing()[2]->mass(),
relativeError()));
}
void StoSFFDecayer::doinit() {
GeneralThreeBodyDecayer::doinit();
if(outgoing().empty()) return;
unsigned int ndiags = getProcessInfo().size();
sca_.resize(ndiags);
fer_.resize(ndiags);
RSfer_.resize(ndiags);
vec_.resize(ndiags);
ten_.resize(ndiags);
four_.resize(ndiags);
for(unsigned int ix = 0;ix < ndiags; ++ix) {
TBDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
// four point vertex
if(!offshell) {
four_[ix] = dynamic_ptr_cast<AbstractFFSSVertexPtr>(current.vertices.first);
continue;
}
if( offshell->CC() ) offshell = offshell->CC();
if(offshell->iSpin() == PDT::Spin0) {
AbstractSSSVertexPtr vert1 = dynamic_ptr_cast<AbstractSSSVertexPtr>
(current.vertices.first);
AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in StoSFFDecayer::doinit()"
<< Exception::runerror;
sca_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1Half) {
AbstractFFSVertexPtr vert1 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.first);
AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a fermion diagram in StoSFFDecayer::doinit()"
<< Exception::runerror;
fer_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractVSSVertexPtr vert1 = dynamic_ptr_cast<AbstractVSSVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a vector diagram in StoSFFDecayer::doinit()"
<< Exception::runerror;
vec_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin2) {
AbstractSSTVertexPtr vert1 = dynamic_ptr_cast<AbstractSSTVertexPtr>
(current.vertices.first);
AbstractFFTVertexPtr vert2 = dynamic_ptr_cast<AbstractFFTVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a tensor diagram in StoSFFDecayer::doinit()"
<< Exception::runerror;
ten_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin3Half) {
AbstractRFSVertexPtr vert1 = dynamic_ptr_cast<AbstractRFSVertexPtr>
(current.vertices.first);
AbstractRFSVertexPtr vert2 = dynamic_ptr_cast<AbstractRFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a RS fermion diagram in StoSFFDecayer::doinit()"
<< Exception::runerror;
RSfer_[ix] = make_pair(vert1, vert2);
}
}
}
-double StoSFFDecayer::me2(const int ichan, const Particle & inpart,
- const ParticleVector & decay,
+void StoSFFDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ ScalarWaveFunction::
+ constructSpinInfo(const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true);
+ for(unsigned int ix=0;ix<decay.size();++ix) {
+ if(decay[ix]->dataPtr()->iSpin()==PDT::Spin0) {
+ ScalarWaveFunction::constructSpinInfo(decay[ix],Helicity::outgoing,true);
+ }
+ else {
+ SpinorWaveFunction::
+ constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
+ }
+ }
+}
+
+double StoSFFDecayer::me2(const int ichan, const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
// particle or CC of particle
- bool cc = (*getProcessInfo().begin()).incoming != inpart.id();
+ bool cc = (*getProcessInfo().begin()).incoming != part.id();
// special handling or first/last call
if(meopt==Initialize) {
ScalarWaveFunction::
- calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming);
- swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),
+ swave_ = ScalarWaveFunction(part.momentum(),part.dataPtr(),
Helicity::incoming);
// fix rho if no correlations
fixRho(rho_);
}
- if(meopt==Terminate) {
- ScalarWaveFunction::
- constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true);
- for(unsigned int ix=0;ix<decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin()==PDT::Spin0) {
- ScalarWaveFunction::constructSpinInfo(decay[ix],Helicity::outgoing,true);
- }
- else {
- SpinorWaveFunction::
- constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
- }
- }
- return 0.;
- }
// get the wavefunctions for all the particles
ScalarWaveFunction outScalar;
unsigned int isca(0);
- for(unsigned int ix=0;ix<decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin()==PDT::Spin0) {
+ for(unsigned int ix=0;ix<outgoing.size();++ix) {
+ if(outgoing[ix]->iSpin()==PDT::Spin0) {
isca = ix;
- outScalar = ScalarWaveFunction(decay[ix]->momentum(),
- decay[ix]->dataPtr(),Helicity::outgoing);
+ outScalar = ScalarWaveFunction(momenta[ix],outgoing[ix],Helicity::outgoing);
}
else {
SpinorWaveFunction::
- calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing);
+ calculateWaveFunctions(outspin_[ix].first,momenta[ix],outgoing[ix],Helicity::outgoing);
outspin_[ix].second.resize(2);
if(outspin_[ix].first[0].wave().Type() == SpinorType::u) {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].first[iy].conjugate();
}
}
else {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].second[iy].conjugate();
}
}
}
}
const vector<vector<double> > cfactors(getColourFactors());
const vector<vector<double> > nfactors(getLargeNcColourFactors());
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
const size_t ncf(numberOfFlows());
vector<Complex> flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.));
vector<GeneralDecayMEPtr>
mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0,
isca==0 ? PDT::Spin0 : PDT::Spin1Half,
isca==1 ? PDT::Spin0 : PDT::Spin1Half,
isca==2 ? PDT::Spin0 : PDT::Spin1Half)));
vector<GeneralDecayMEPtr>
mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin0,
isca == 0 ? PDT::Spin0 : PDT::Spin1Half,
isca == 1 ? PDT::Spin0 : PDT::Spin1Half,
isca == 2 ? PDT::Spin0 : PDT::Spin1Half)));
static const unsigned int out2[3]={1,0,0},out3[3]={2,2,1};
for(unsigned int s1 = 0;s1 < 2; ++s1) {
for(unsigned int s2 = 0;s2 < 2; ++s2) {
flows = vector<Complex>(ncf, Complex(0.));
largeflows = vector<Complex>(ncf, Complex(0.));
unsigned int idiag(0);
for(vector<TBDiagram>::const_iterator dit = getProcessInfo().begin();
dit != getProcessInfo().end(); ++dit) {
// channels if selecting
if( ichan >= 0 && diagramMap()[ichan] != idiag ) {
++idiag;
continue;
}
tcPDPtr offshell = dit->intermediate;
Complex diag;
if(offshell) {
if(cc&&offshell->CC()) offshell=offshell->CC();
double sign = out3[dit->channelType] < out2[dit->channelType] ? 1. : -1.;
// intermediate scalar
if (offshell->iSpin() == PDT::Spin0) {
ScalarWaveFunction inters = sca_[idiag].first->
evaluate(scale, widthOption(), offshell, swave_, outScalar);
unsigned int h1(s1),h2(s2);
if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2);
- if(decay[out2[dit->channelType]]->id()<0&&
- decay[out3[dit->channelType]]->id()>0) {
+ if(outgoing[out2[dit->channelType]]->id()<0&&
+ outgoing[out3[dit->channelType]]->id()>0) {
diag =-sign*sca_[idiag].second->
evaluate(scale,
outspin_[out2[dit->channelType]].first [h1],
outspin_[out3[dit->channelType]].second[h2],inters);
}
else {
diag = sign*sca_[idiag].second->
evaluate(scale,
outspin_[out3[dit->channelType]].first [h2],
outspin_[out2[dit->channelType]].second[h1],inters);
}
}
// intermediate fermion
else if(offshell->iSpin() == PDT::Spin1Half) {
int iferm =
- decay[out2[dit->channelType]]->dataPtr()->iSpin()==PDT::Spin1Half
+ outgoing[out2[dit->channelType]]->iSpin()==PDT::Spin1Half
? out2[dit->channelType] : out3[dit->channelType];
unsigned int h1(s1),h2(s2);
if(dit->channelType>iferm) swap(h1,h2);
sign = iferm<dit->channelType ? 1. : -1.;
- if((decay[dit->channelType]->id() < 0 &&decay[iferm]->id() > 0 ) ||
- (decay[dit->channelType]->id()*offshell->id()>0)) {
+ if((outgoing[dit->channelType]->id() < 0 &&outgoing[iferm]->id() > 0 ) ||
+ (outgoing[dit->channelType]->id()*offshell->id()>0)) {
SpinorWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].first [h1],swave_);
diag = -sign*fer_[idiag].second->
evaluate(scale,inters,outspin_[iferm].second[h2],outScalar);
}
else {
SpinorBarWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].second[h1],swave_);
diag = sign*fer_[idiag].second->
evaluate(scale,outspin_[iferm].first [h2],inters,outScalar);
}
}
// intermediate vector
else if(offshell->iSpin() == PDT::Spin1) {
VectorWaveFunction interv = vec_[idiag].first->
evaluate(scale, widthOption(), offshell, swave_, outScalar);
unsigned int h1(s1),h2(s2);
if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2);
- if(decay[out2[dit->channelType]]->id()<0&&
- decay[out3[dit->channelType]]->id()>0) {
+ if(outgoing[out2[dit->channelType]]->id()<0&&
+ outgoing[out3[dit->channelType]]->id()>0) {
diag =-sign*vec_[idiag].second->
evaluate(scale,
outspin_[out2[dit->channelType]].first [h1],
outspin_[out3[dit->channelType]].second[h2],interv);
}
else {
diag = sign*vec_[idiag].second->
evaluate(scale,
outspin_[out3[dit->channelType]].first [h2],
outspin_[out2[dit->channelType]].second[h1],interv);
}
}
// intermediate tensor
else if(offshell->iSpin() == PDT::Spin2) {
TensorWaveFunction intert = ten_[idiag].first->
evaluate(scale, widthOption(), offshell, swave_, outScalar);
unsigned int h1(s1),h2(s2);
if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2);
- if(decay[out2[dit->channelType]]->id()<0&&
- decay[out3[dit->channelType]]->id()>0) {
+ if(outgoing[out2[dit->channelType]]->id()<0&&
+ outgoing[out3[dit->channelType]]->id()>0) {
diag =-sign*ten_[idiag].second->
evaluate(scale,
outspin_[out2[dit->channelType]].first [h1],
outspin_[out3[dit->channelType]].second[h2],intert);
}
else {
diag = sign*ten_[idiag].second->
evaluate(scale,
outspin_[out3[dit->channelType]].first [h2],
outspin_[out2[dit->channelType]].second[h1],intert);
}
}
// intermediate RS fermion
else if(offshell->iSpin() == PDT::Spin3Half) {
int iferm =
- decay[out2[dit->channelType]]->dataPtr()->iSpin()==PDT::Spin1Half
+ outgoing[out2[dit->channelType]]->iSpin()==PDT::Spin1Half
? out2[dit->channelType] : out3[dit->channelType];
unsigned int h1(s1),h2(s2);
if(dit->channelType>iferm) swap(h1,h2);
sign = iferm<dit->channelType ? 1. : -1.;
- if((decay[dit->channelType]->id() < 0 &&decay[iferm]->id() > 0 ) ||
- (decay[dit->channelType]->id()*offshell->id()>0)) {
+ if((outgoing[dit->channelType]->id() < 0 &&outgoing[iferm]->id() > 0 ) ||
+ (outgoing[dit->channelType]->id()*offshell->id()>0)) {
RSSpinorWaveFunction inters = RSfer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].first [h1],swave_);
diag = -sign*RSfer_[idiag].second->
evaluate(scale,inters,outspin_[iferm].second[h2],outScalar);
}
else {
RSSpinorBarWaveFunction inters = RSfer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].second[h1],swave_);
diag = sign*RSfer_[idiag].second->
evaluate(scale,outspin_[iferm].first [h2],inters,outScalar);
}
}
// unknown
else throw Exception()
<< "Unknown intermediate in StoSFFDecayer::me2()"
<< Exception::runerror;
}
// four point diagram
else {
unsigned int o2 = isca > 0 ? 0 : 1;
unsigned int o3 = isca < 2 ? 2 : 1;
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag =-four_[idiag]->
evaluate(scale, outspin_[o2].first[s1],
outspin_[o3].second[s2], outScalar, swave_);
}
else {
diag = four_[idiag]->
evaluate(scale, outspin_[o3].first[s2],
outspin_[o2].second[s1], outScalar, swave_);
}
}
// matrix element for the different colour flows
if(ichan < 0) {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
else {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1 != colourFlow()) continue;
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
++idiag;
}
for(unsigned int ix = 0; ix < ncf; ++ix) {
if(isca == 0) {
(*mes[ix])(0, 0, s1, s2) = flows[ix];
(*mel[ix])(0, 0, s1, s2) = largeflows[ix];
}
else if(isca == 1 ) {
(*mes[ix])(0, s1, 0, s2) = flows[ix];
(*mel[ix])(0, s1, 0, s2) = largeflows[ix];
}
else if(isca == 2) {
(*mes[ix])(0, s1,s2, 0) = flows[ix];
(*mel[ix])(0, s1,s2, 0) = largeflows[ix] ;
}
}
}
}
double me2(0.);
if(ichan < 0) {
vector<double> pflows(ncf,0.);
for(unsigned int ix = 0; ix < ncf; ++ix) {
for(unsigned int iy = 0; iy < ncf; ++ iy) {
double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real();
me2 += con;
if(ix == iy) {
con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real();
pflows[ix] += con;
}
}
}
double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.));
ptotal *= UseRandom::rnd();
for(unsigned int ix = 0;ix < pflows.size(); ++ix) {
if(ptotal <= pflows[ix]) {
colourFlow(ix);
ME(mes[ix]);
break;
}
ptotal -= pflows[ix];
}
}
else {
unsigned int iflow = colourFlow();
me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real();
}
// return the matrix element squared
return me2;
}
diff --git a/Decay/General/StoSFFDecayer.h b/Decay/General/StoSFFDecayer.h
--- a/Decay/General/StoSFFDecayer.h
+++ b/Decay/General/StoSFFDecayer.h
@@ -1,163 +1,172 @@
// -*- C++ -*-
#ifndef THEPEG_StoSFFDecayer_H
#define THEPEG_StoSFFDecayer_H
//
// This is the declaration of the StoSFFDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSSVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* The StoSFFDecayer class provides the general matrix element for
* scalar decays into another scalar and a fermion-antifermion pair.
*
* @see \ref StoSFFDecayerInterfaces "The interfaces"
* defined for StoSFFDecayer.
*/
class StoSFFDecayer: public GeneralThreeBodyDecayer {
public:
-
+
/**
- * 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;
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
StoSFFDecayer & operator=(const StoSFFDecayer &);
private:
/**
* Store the vertices for scalar intermediate
*/
vector<pair<AbstractSSSVertexPtr, AbstractFFSVertexPtr> > sca_;
/**
* Store the vertices for spin-\f$\frac12\f$ fermion intermediate
*/
vector<pair<AbstractFFSVertexPtr, AbstractFFSVertexPtr> > fer_;
/**
* Store the vertices for spin-\f$\frac32\f$ fermion intermediate
*/
vector<pair<AbstractRFSVertexPtr, AbstractRFSVertexPtr> > RSfer_;
/**
* Store the vertices for vector intermediate
*/
vector<pair<AbstractVSSVertexPtr, AbstractFFVVertexPtr> > vec_;
/**
* Store the vertices for tensor intermediate
*/
vector<pair<AbstractSSTVertexPtr, AbstractFFTVertexPtr> > ten_;
/**
* Store the vertices for four point diagrams
*/
vector<AbstractFFSSVertexPtr> four_;
/**
* Spin density matrix
*/
mutable RhoDMatrix rho_;
/**
* Scalar wavefunction
*/
mutable ScalarWaveFunction swave_;
/**
* Spinor wavefunctions
*/
mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > outspin_[3];
};
}
#endif /* THEPEG_StoSFFDecayer_H */
diff --git a/Decay/General/VtoFFVDecayer.cc b/Decay/General/VtoFFVDecayer.cc
--- a/Decay/General/VtoFFVDecayer.cc
+++ b/Decay/General/VtoFFVDecayer.cc
@@ -1,371 +1,372 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the VtoFFVDecayer class.
//
#include "VtoFFVDecayer.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/PDT/ThreeBodyAllOnCalculator.h"
#include "Herwig/Decay/GeneralDecayMatrixElement.h"
#include <numeric>
using namespace Herwig;
using namespace ThePEG;
using namespace ThePEG::Helicity;
IBPtr VtoFFVDecayer::clone() const {
return new_ptr(*this);
}
IBPtr VtoFFVDecayer::fullclone() const {
return new_ptr(*this);
}
void VtoFFVDecayer::persistentOutput(PersistentOStream & os) const {
os << sca_ << fer_ << vec_ << ten_;
}
void VtoFFVDecayer::persistentInput(PersistentIStream & is, int) {
is >> sca_ >> fer_ >> vec_ >> ten_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<VtoFFVDecayer,GeneralThreeBodyDecayer>
describeHerwigVtoFFVDecayer("Herwig::VtoFFVDecayer", "Herwig.so");
void VtoFFVDecayer::Init() {
static ClassDocumentation<VtoFFVDecayer> documentation
("The VtoFFVDecayer class implements the general three-body "
"decay of a vector to a two fermions and a vector.");
}
WidthCalculatorBasePtr VtoFFVDecayer::
threeBodyMEIntegrator(const DecayMode & ) const {
vector<int> intype;
vector<Energy> inmass,inwidth;
vector<double> inpow,inweights;
constructIntegratorChannels(intype,inmass,inwidth,inpow,inweights);
return new_ptr(ThreeBodyAllOnCalculator<VtoFFVDecayer>
(inweights,intype,inmass,inwidth,inpow,*this,0,
outgoing()[0]->mass(),outgoing()[1]->mass(),
outgoing()[2]->mass(),relativeError()));
}
void VtoFFVDecayer::doinit() {
GeneralThreeBodyDecayer::doinit();
if(outgoing().empty()) return;
unsigned int ndiags = getProcessInfo().size();
sca_.resize(ndiags);
fer_.resize(ndiags);
vec_.resize(ndiags);
ten_.resize(ndiags);
for(unsigned int ix = 0;ix < ndiags; ++ix) {
TBDiagram current = getProcessInfo()[ix];
tcPDPtr offshell = current.intermediate;
if( offshell->CC() ) offshell = offshell->CC();
if(offshell->iSpin() == PDT::Spin0) {
AbstractVVSVertexPtr vert1 = dynamic_ptr_cast<AbstractVVSVertexPtr>
(current.vertices.first);
AbstractFFSVertexPtr vert2 = dynamic_ptr_cast<AbstractFFSVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a scalar diagram in VtoFFVDecayer::doinit()"
<< Exception::runerror;
sca_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1Half) {
AbstractFFVVertexPtr vert1 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a fermion diagram in VtoFFVDecayer::doinit()"
<< Exception::runerror;
fer_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin1) {
AbstractVVVVertexPtr vert1 = dynamic_ptr_cast<AbstractVVVVertexPtr>
(current.vertices.first);
AbstractFFVVertexPtr vert2 = dynamic_ptr_cast<AbstractFFVVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a vector diagram in VtoFFVDecayer::doinit()"
<< Exception::runerror;
vec_[ix] = make_pair(vert1, vert2);
}
else if(offshell->iSpin() == PDT::Spin2) {
AbstractVVTVertexPtr vert1 = dynamic_ptr_cast<AbstractVVTVertexPtr>
(current.vertices.first);
AbstractFFTVertexPtr vert2 = dynamic_ptr_cast<AbstractFFTVertexPtr>
(current.vertices.second);
if(!vert1||!vert2) throw Exception()
<< "Invalid vertices for a tensor diagram in VtoFFVDecayer::doinit()"
<< Exception::runerror;
ten_[ix] = make_pair(vert1, vert2);
}
}
}
+void VtoFFVDecayer::
+constructSpinInfo(const Particle & part, ParticleVector decay) const {
+ VectorWaveFunction::
+ constructSpinInfo(inVector_,const_ptr_cast<tPPtr>(&part),
+ Helicity::incoming,true,false);
+ for(unsigned int ix=0;ix<decay.size();++ix) {
+ if(decay[ix]->dataPtr()->iSpin()==PDT::Spin1) {
+ VectorWaveFunction::constructSpinInfo(outVector_,decay[ix],
+ Helicity::outgoing,true,false);
+ }
+ else {
+ SpinorWaveFunction::
+ constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
+ }
+ }
+}
-double VtoFFVDecayer::me2(const int ichan, const Particle & inpart,
- const ParticleVector & decay,
+double VtoFFVDecayer::me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
MEOption meopt) const {
// particle or CC of particle
- bool cc = (*getProcessInfo().begin()).incoming != inpart.id();
+ bool cc = (*getProcessInfo().begin()).incoming != part.id();
// special handling or first/last call
if(meopt==Initialize) {
VectorWaveFunction::
- calculateWaveFunctions(inVector_,rho_,const_ptr_cast<tPPtr>(&inpart),
+ calculateWaveFunctions(inVector_,rho_,const_ptr_cast<tPPtr>(&part),
Helicity::incoming,false);
// fix rho if no correlations
fixRho(rho_);
}
-
- if(meopt==Terminate) {
- VectorWaveFunction::
- constructSpinInfo(inVector_,const_ptr_cast<tPPtr>(&inpart),
- Helicity::incoming,true,false);
- for(unsigned int ix=0;ix<decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin()==PDT::Spin1) {
- VectorWaveFunction::constructSpinInfo(outVector_,decay[ix],
- Helicity::outgoing,true,false);
- }
- else {
- SpinorWaveFunction::
- constructSpinInfo(outspin_[ix].first,decay[ix],Helicity::outgoing,true);
- }
- }
- }
unsigned int ivec(0);
bool massless = false;
- for(unsigned int ix = 0; ix < decay.size();++ix) {
- if(decay[ix]->dataPtr()->iSpin() == PDT::Spin1) {
- massless = decay[ix]->id()==ParticleID::g || decay[ix]->id()==ParticleID::gamma;
+ for(unsigned int ix = 0; ix < outgoing.size();++ix) {
+ if(outgoing[ix]->iSpin() == PDT::Spin1) {
+ massless = outgoing[ix]->id()==ParticleID::g || outgoing[ix]->id()==ParticleID::gamma;
ivec = ix;
VectorWaveFunction::
- calculateWaveFunctions(outVector_, decay[ix], Helicity::outgoing, massless );
+ calculateWaveFunctions(outVector_, momenta[ix], outgoing[ix], Helicity::outgoing, massless );
}
else {
SpinorWaveFunction::
- calculateWaveFunctions(outspin_[ix].first,decay[ix],Helicity::outgoing);
+ calculateWaveFunctions(outspin_[ix].first,momenta[ix], outgoing[ix],Helicity::outgoing);
outspin_[ix].second.resize(2);
// Need a ubar and a v spinor
if(outspin_[ix].first[0].wave().Type() == SpinorType::u) {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].first[iy].conjugate();
}
}
else {
for(unsigned int iy = 0; iy < 2; ++iy) {
outspin_[ix].second[iy] = outspin_[ix].first[iy].bar();
outspin_[ix].second[iy].conjugate();
}
}
}
}
const vector<vector<double> > cfactors(getColourFactors());
const vector<vector<double> > nfactors(getLargeNcColourFactors());
- Energy2 scale(sqr(inpart.mass()));
+ Energy2 scale(sqr(part.mass()));
const size_t ncf(numberOfFlows());
vector<Complex> flows(ncf, Complex(0.)), largeflows(ncf, Complex(0.));
// setup the DecayMatrixElement
vector<GeneralDecayMEPtr>
mes(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1,
ivec == 0 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 1 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 2 ? PDT::Spin1 : PDT::Spin1Half)));
vector<GeneralDecayMEPtr>
mel(ncf,new_ptr(GeneralDecayMatrixElement(PDT::Spin1,
ivec == 0 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 1 ? PDT::Spin1 : PDT::Spin1Half,
ivec == 2 ? PDT::Spin1 : PDT::Spin1Half)));
//the channel possiblities
static const unsigned int out2[3] = {1,0,0}, out3[3] = {2,2,1};
for(unsigned int vi = 0; vi < 3; ++vi) {
for(unsigned int s1 = 0; s1 < 2; ++s1) {
for(unsigned int s2 = 0; s2 < 2; ++s2) {
for(unsigned int v1 = 0; v1 < 3; ++v1) {
if ( massless && v1 == 1 ) continue;
flows = vector<Complex>(ncf, Complex(0.));
largeflows = vector<Complex>(ncf, Complex(0.));
unsigned int idiag(0);
for(vector<TBDiagram>::const_iterator dit=getProcessInfo().begin();
dit!=getProcessInfo().end();++dit) {
// channels if selecting
if( ichan >=0 && diagramMap()[ichan] != idiag ) {
++idiag;
continue;
}
tcPDPtr offshell = dit->intermediate;
if(cc&&offshell->CC()) offshell=offshell->CC();
Complex diag;
unsigned int o2(out2[dit->channelType]), o3(out3[dit->channelType]);
double sign = (o3 < o2) ? 1. : -1.;
// intermediate scalar
if(offshell->iSpin() == PDT::Spin0) {
ScalarWaveFunction inters = sca_[idiag].first->
evaluate(scale, widthOption(), offshell, outVector_[v1],
inVector_[vi]);
unsigned int h1(s1),h2(s2);
if(o2 > o3) swap(h1, h2);
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag = -sign*sca_[idiag].second->
evaluate(scale,outspin_[o2].first[h1],
outspin_[o3].second[h2],inters);
}
else {
diag = sign*sca_[idiag].second->
evaluate(scale, outspin_[o3].first [h2],
outspin_[o2].second[h1],inters);
}
}
// intermediate fermion
else if(offshell->iSpin() == PDT::Spin1Half) {
- int iferm = (decay[o2]->dataPtr()->iSpin() == PDT::Spin1Half)
+ int iferm = (outgoing[o2]->iSpin() == PDT::Spin1Half)
? o2 : o3;
unsigned int h1(s1),h2(s2);
if(dit->channelType > iferm) swap(h1, h2);
sign = iferm < dit->channelType ? 1. : -1.;
- if(decay[dit->channelType]->id() < 0 && decay[iferm]->id() > 0 ) {
+ if(outgoing[dit->channelType]->id() < 0 && outgoing[iferm]->id() > 0 ) {
SpinorWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].first[h1], inVector_[vi]);
diag = -sign*fer_[idiag].second->
evaluate(scale,inters,outspin_[iferm].second[h2], outVector_[v1]);
}
else {
SpinorBarWaveFunction inters = fer_[idiag].first->
evaluate(scale,widthOption(),offshell,
outspin_[dit->channelType].second[h1],inVector_[vi]);
diag = sign*fer_[idiag].second->
evaluate(scale,outspin_[iferm].first [h2],inters, outVector_[v1]);
}
}
// intermediate vector
else if(offshell->iSpin() == PDT::Spin1) {
VectorWaveFunction interv = vec_[idiag].first->
evaluate(scale, widthOption(), offshell, outVector_[v1],
inVector_[vi]);
unsigned int h1(s1),h2(s2);
if(o2 > o3) swap(h1,h2);
- if(decay[o2]->id() < 0 && decay[o3]->id() > 0) {
+ if(outgoing[o2]->id() < 0 && outgoing[o3]->id() > 0) {
diag =-sign*vec_[idiag].second->
evaluate(scale, outspin_[o2].first[h1],
outspin_[o3].second[h2], interv);
}
else {
diag = sign*vec_[idiag].second->
evaluate(scale, outspin_[o3].first[h2],
outspin_[o2].second[h1], interv);
}
}
else if(offshell->iSpin() == PDT::Spin2) {
TensorWaveFunction intert = ten_[idiag].first->
evaluate(scale, widthOption(), offshell, inVector_[vi],
outVector_[v1]);
unsigned int h1(s1),h2(s2);
if(out2[dit->channelType]>out3[dit->channelType]) swap(h1,h2);
- if(decay[out2[dit->channelType]]->id()<0&&
- decay[out3[dit->channelType]]->id()>0) {
+ if(outgoing[out2[dit->channelType]]->id()<0&&
+ outgoing[out3[dit->channelType]]->id()>0) {
diag =-sign*ten_[idiag].second->
evaluate(scale,
outspin_[out2[dit->channelType]].first [h1],
outspin_[out3[dit->channelType]].second[h2],intert);
}
else {
diag = sign*ten_[idiag].second->
evaluate(scale,
outspin_[out3[dit->channelType]].first [h2],
outspin_[out2[dit->channelType]].second[h1],intert);
}
}
// unknown
else throw Exception()
<< "Unknown intermediate in VtoFFVDecayer::me2()"
<< Exception::runerror;
// matrix element for the different colour flows
if(ichan < 0) {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
else {
for(unsigned iy = 0; iy < dit->colourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1 != colourFlow()) continue;
flows[dit->colourFlow[iy].first - 1] +=
dit->colourFlow[iy].second * diag;
}
for(unsigned iy = 0; iy < dit->largeNcColourFlow.size(); ++iy) {
if(dit->colourFlow[iy].first - 1!=colourFlow()) continue;
largeflows[dit->largeNcColourFlow[iy].first - 1] +=
dit->largeNcColourFlow[iy].second * diag;
}
}
++idiag;
} //end of diagrams
// now add the flows to the me2 with appropriate colour factors
for(unsigned int ix = 0; ix < ncf; ++ix) {
if (ivec == 0) {
(*mes[ix])(vi, v1, s1, s2) = flows[ix];
(*mel[ix])(vi, v1, s1, s2) = largeflows[ix];
}
else if(ivec == 1) {
(*mes[ix])(vi, s1, v1, s2) = flows[ix];
(*mel[ix])(vi, s1, v1, s2) = largeflows[ix];
}
else if(ivec == 2) {
(*mes[ix])(vi, s1, s2, v1) = flows[ix];
(*mel[ix])(vi, s1, s2, v1) = largeflows[ix];
}
}
}
}
}
}
double me2(0.);
if(ichan < 0) {
vector<double> pflows(ncf,0.);
for(unsigned int ix = 0; ix < ncf; ++ix) {
for(unsigned int iy = 0; iy < ncf; ++ iy) {
double con = cfactors[ix][iy]*(mes[ix]->contract(*mes[iy],rho_)).real();
me2 += con;
if(ix == iy) {
con = nfactors[ix][iy]*(mel[ix]->contract(*mel[iy],rho_)).real();
pflows[ix] += con;
}
}
}
double ptotal(std::accumulate(pflows.begin(),pflows.end(),0.));
ptotal *= UseRandom::rnd();
for(unsigned int ix = 0;ix < pflows.size(); ++ix) {
if(ptotal <= pflows[ix]) {
colourFlow(ix);
ME(mes[ix]);
break;
}
ptotal -= pflows[ix];
}
}
else {
unsigned int iflow = colourFlow();
me2 = nfactors[iflow][iflow]*(mel[iflow]->contract(*mel[iflow],rho_)).real();
}
// return the matrix element squared
return me2;
}
diff --git a/Decay/General/VtoFFVDecayer.h b/Decay/General/VtoFFVDecayer.h
--- a/Decay/General/VtoFFVDecayer.h
+++ b/Decay/General/VtoFFVDecayer.h
@@ -1,161 +1,170 @@
// -*- C++ -*-
#ifndef THEPEG_VtoFFVDecayer_H
#define THEPEG_VtoFFVDecayer_H
//
// This is the declaration of the VtoFFVDecayer class.
//
#include "GeneralThreeBodyDecayer.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the VtoFFVDecayer class.
*
* @see \ref VtoFFVDecayerInterfaces "The interfaces"
* defined for VtoFFVDecayer.
*/
class VtoFFVDecayer: public GeneralThreeBodyDecayer {
public:
+
+ /**
+ * Return the matrix element squared for a given mode and phase-space channel.
+ * @param ichan The channel we are calculating the matrix element for.
+ * @param part The decaying Particle.
+ * @param outgoing The particles produced in the decay
+ * @param momenta The momenta of the particles produced in the decay
+ * @param meopt Option for the calculation of the matrix element
+ * @return The matrix element squared for the phase-space configuration.
+ */
+ double me2(const int ichan,const Particle & part,
+ const tPDVector & outgoing,
+ const vector<Lorentz5Momentum> & momenta,
+ MEOption meopt) const;
/**
- * 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
- * @return The matrix element squared for the phase-space configuration.
+ * Construct the SpinInfos for the particles produced in the decay
*/
- virtual double me2(const int ichan, const Particle & part,
- const ParticleVector & decay, MEOption meopt) const;
+ virtual void constructSpinInfo(const Particle & part,
+ ParticleVector outgoing) const;
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of
* calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
public:
/** @name Functions used by the persistent I/O system. */
//@{
/**
* Function used to write out object persistently.
* @param os the persistent output stream written to.
*/
void persistentOutput(PersistentOStream & os) const;
/**
* Function used to read in object persistently.
* @param is the persistent input stream read from.
* @param version the version number of the object when written.
*/
void persistentInput(PersistentIStream & is, int version);
//@}
/**
* The standard Init function used to initialize the interfaces.
* Called exactly once for each class by the class description system
* before the main function starts or
* when this class is dynamically loaded.
*/
static void Init();
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
VtoFFVDecayer & operator=(const VtoFFVDecayer &);
private:
/**
* Store the vertices for scalar intrermediate
*/
vector<pair<AbstractVVSVertexPtr, AbstractFFSVertexPtr> > sca_;
/**
* Store the vertices for fermion intrermediate
*/
vector<pair<AbstractFFVVertexPtr, AbstractFFVVertexPtr> > fer_;
/**
* Store the vertices for vector intrermediate
*/
vector<pair<AbstractVVVVertexPtr, AbstractFFVVertexPtr> > vec_;
/**
* Store the vertices for vector intrermediate
*/
vector<pair<AbstractVVTVertexPtr, AbstractFFTVertexPtr> > ten_;
/**
* Spinr density matrix
*/
mutable RhoDMatrix rho_;
/**
* Polarization vectors for the decaying particle
*/
mutable vector<VectorWaveFunction> inVector_;
/**
* Scalar wavefunction for the decay products
*/
mutable ScalarWaveFunction swave_;
/**
* Polarization vectors for the decay products
*/
mutable vector<VectorWaveFunction> outVector_;
/**
* Spinors for the decay products
*/
mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > outspin_[3];
};
}
#endif /* THEPEG_VtoFFVDecayer_H */
diff --git a/Decay/PhaseSpaceMode.cc b/Decay/PhaseSpaceMode.cc
--- a/Decay/PhaseSpaceMode.cc
+++ b/Decay/PhaseSpaceMode.cc
@@ -1,491 +1,492 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the PhaseSpaceMode class.
//
#include "PhaseSpaceMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/EventRecord/Particle.h"
#include "ThePEG/Repository/UseRandom.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig/Utilities/Kinematics.h"
using namespace Herwig;
using namespace ThePEG::Helicity;
void PhaseSpaceMode::persistentOutput(PersistentOStream & os) const {
os << incoming_ << channels_ << maxWeight_ << outgoing_ << outgoingCC_
<< partial_ << widthGen_ << massGen_ << testOnShell_
<< ounit(eMax_,GeV) << nRand_;
}
void PhaseSpaceMode::persistentInput(PersistentIStream & is, int) {
is >> incoming_ >> channels_ >> maxWeight_ >> outgoing_ >> outgoingCC_
>> partial_ >> widthGen_ >> massGen_ >> testOnShell_
>> iunit(eMax_,GeV) >> nRand_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<PhaseSpaceMode,Base>
describeHerwigPhaseSpaceMode("Herwig::PhaseSpaceMode", "PhaseSpaceMode.so");
void PhaseSpaceMode::Init() {
static ClassDocumentation<PhaseSpaceMode> documentation
("The PhaseSpaceMode classs contains a number of phase space"
" channels for the integration of a particular decay mode");
}
ParticleVector
PhaseSpaceMode::generateDecay(const Particle & inpart,
tcDecayIntegrator2Ptr decayer,
bool intermediates,bool cc) {
decayer->ME(DecayMEPtr());
eps_=decayer->eps_;
// compute the prefactor
Energy prewid = (widthGen_&&partial_>=0) ?
widthGen_->partialWidth(partial_,inpart.mass()) :
incoming_.first->width();
InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
// boosts to/from rest
Boost bv =-inpart.momentum().boostVector();
double gammarest = inpart.momentum().e()/inpart.momentum().mass();
LorentzRotation boostToRest( bv,gammarest);
LorentzRotation boostFromRest(-bv,gammarest);
// construct a new particle which is at rest
Particle inrest(inpart);
inrest.transform(boostToRest);
double wgt(0.);
vector<Lorentz5Momentum> momenta(outgoing_.size());
int ichan;
unsigned int ncount(0);
try {
do {
// phase-space pieces of the weight
fillStack();
wgt = pre*weight(ichan,inrest.momentum(),momenta);
// matrix element piece
wgt *= decayer->me2(-1,inrest,!cc ? outgoing_ : outgoingCC_,
momenta,
ncount ==0 ? DecayIntegrator2::Initialize : DecayIntegrator2::Calculate);
++ncount;
if(wgt>maxWeight_) {
CurrentGenerator::log() << "Resetting max weight for decay "
<< inrest.PDGName() << " -> ";
for(tcPDPtr part : outgoing_)
CurrentGenerator::log() << " " << part->PDGName();
CurrentGenerator::log() << " " << maxWeight_ << " " << wgt
<< " " << inrest.mass()/MeV << "\n";
maxWeight_=wgt;
}
}
while(maxWeight_*UseRandom::rnd()>wgt&&ncount<decayer->nTry_);
}
catch (Veto) {
// restore the incoming particle to its original state
inrest.transform(boostFromRest);
throw Veto();
}
if(ncount>=decayer->nTry_) {
CurrentGenerator::log() << "The decay " << inrest.PDGName() << " -> ";
for(tcPDPtr part : outgoing_)
CurrentGenerator::log() << " " << part->PDGName();
CurrentGenerator::log() << " " << maxWeight_ << " " << decayer->nTry_
<< " is too inefficient for the particle "
<< inpart << "vetoing the decay \n";
momenta.clear();
throw Veto();
}
// construct the particles
ParticleVector output(outgoing_.size());
for(unsigned int ix=0;ix<outgoing_.size();++ix)
output[ix] = (!cc ? outgoing_[ix] : outgoingCC_[ix] )->produceParticle(momenta[ix]);
// set up the spin correlations
decayer->constructSpinInfo(inrest,output);
const_ptr_cast<tPPtr>(&inpart)->spinInfo(inrest.spinInfo());
constructVertex(inpart,output,decayer);
// return if intermediate particles not required
if(channels_.empty()||!intermediates) {
for(tPPtr part : output) part->transform(boostFromRest);
}
// find the intermediate particles
else {
// select the channel
vector<double> mewgts(channels_.size(),0.0);
double total=0.;
for(unsigned int ix=0,N=channels_.size();ix<N;++ix) {
mewgts[ix]=decayer->me2(ix,inrest,!cc ? outgoing_ : outgoingCC_,
momenta,DecayIntegrator2::Calculate);
total+=mewgts[ix];
}
// randomly pick a channel
total *= UseRandom::rnd();
int iChannel = -1;
do {
++iChannel;
total-=mewgts[iChannel];
}
while(iChannel<int(channels_.size()) && total>0.);
+ iChannel_ = iChannel;
// apply boost
for(tPPtr part : output) part->transform(boostFromRest);
channels_[iChannel].generateIntermediates(cc,inpart,output);
}
decayer->ME(DecayMEPtr());
// return particles;
return output;
}
// flat phase space generation and weight
Energy PhaseSpaceMode::flatPhaseSpace(const Lorentz5Momentum & in,
vector<Lorentz5Momentum> & momenta,
bool onShell) const {
double wgt(1.);
// masses of the particles
vector<Energy> mass = externalMasses(in.mass(),wgt,onShell);
// two body decay
double ctheta = 2.*rStack_.top() - 1.;
rStack_.pop();
double phi = Constants::twopi*rStack_.top();
rStack_.pop();
if(! Kinematics::twoBodyDecay(in, mass[0], mass[1],
ctheta, phi, momenta[0],momenta[1]))
throw Exception() << "Incoming mass - Outgoing mass negative in "
<< "PhaseSpaceMode::flatPhaseSpace()"
<< Exception::eventerror;
wgt *= Kinematics::pstarTwoBodyDecay(in.mass(),mass[0],mass[1])
/8./Constants::pi/in.mass();
return wgt*in.mass();
}
// generate the masses of the external particles
vector<Energy> PhaseSpaceMode::externalMasses(Energy inmass,double & wgt,
bool onShell) const {
vector<Energy> mass(outgoing_.size());
vector<int> notdone;
Energy mlow(ZERO);
// set masses of stable particles and limits
for(unsigned int ix=0;ix<outgoing_.size();++ix) {
// get the mass of the particle if can't use weight
if(onShell)
mass[ix] = outgoing_[ix]->mass();
else if(!massGen_[ix] || outgoing_[ix]->stable()) {
mass[ix] = outgoing_[ix]->generateMass();
mlow += mass[ix];
}
else {
mass[ix] = ZERO;
notdone.push_back(ix);
mlow+=max(outgoing_[ix]->mass()-outgoing_[ix]->widthLoCut(),ZERO);
}
}
if(mlow>inmass) throw Veto();
// now we need to generate the masses for the particles we haven't
for( ;!notdone.empty();) {
unsigned int iloc=long(UseRandom::rnd()*(notdone.size()-1));
Energy low = max(outgoing_[notdone[iloc]]->mass()-outgoing_[notdone[iloc]]->widthLoCut(),ZERO);
mlow-=low;
double wgttemp;
mass[notdone[iloc]]= massGen_[notdone[iloc]]->mass(wgttemp,*outgoing_[notdone[iloc]],low,inmass-mlow,rStack_.top());
rStack_.pop();
assert(mass[notdone[iloc]]>=low&&mass[notdone[iloc]]<=inmass-mlow);
wgt *= wgttemp;
mlow += mass[notdone[iloc]];
notdone.erase(notdone.begin()+iloc);
}
return mass;
}
// construct the vertex for spin corrections
void PhaseSpaceMode::constructVertex(const Particle & inpart,
const ParticleVector & decay,
tcDecayIntegrator2Ptr decayer) const {
// construct the decay vertex
VertexPtr vertex(new_ptr(DecayVertex()));
DVertexPtr Dvertex(dynamic_ptr_cast<DVertexPtr>(vertex));
// set the incoming particle for the decay vertex
(inpart.spinInfo())->decayVertex(vertex);
for(unsigned int ix=0;ix<decay.size();++ix) {
decay[ix]->spinInfo()->productionVertex(vertex);
}
// set the matrix element
Dvertex->ME(decayer->ME());
decayer->ME(DecayMEPtr());
}
// initialise the phase space
Energy PhaseSpaceMode::initializePhaseSpace(bool init, tcDecayIntegrator2Ptr decayer,
bool onShell) {
decayer->ME(DecayMEPtr());
eps_=decayer->eps_;
Energy output(ZERO);
// ensure that the weights add up to one
if(!channels_.empty()) {
double temp=0.;
for(const PhaseSpaceChannel & channel : channels_) temp+= channel.weight();
for(PhaseSpaceChannel & channel : channels_) channel.weight(channel.weight()/temp);
}
if(!init) return ZERO;
ThePEG::PPtr inpart=incoming_.first->produceParticle();
// now if using flat phase space
maxWeight_=0.;
if(channels_.empty()) {
vector<Lorentz5Momentum> momenta(outgoing_.size());
double wsum=0.,wsqsum=0.;
Energy m0,mmin(ZERO);
for(tcPDPtr part : outgoing_) mmin += part->massMin();
for(unsigned int ix=0;ix<decayer->nPoint_;++ix) {
// set the mass of the decaying particle
m0 = !onShell ? inpart->dataPtr()->generateMass() : inpart->dataPtr()->mass();
double wgt=0.;
if(m0<=mmin) continue;
inpart->set5Momentum(Lorentz5Momentum(m0));
// compute the prefactor
Energy prewid = (widthGen_&&partial_>=0) ?
widthGen_->partialWidth(partial_,inpart->mass()) :
incoming_.first->width();
InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
// generate the weight for this point
try {
int dummy;
// phase-space piece
fillStack();
wgt = pre*weight(dummy,inpart->momentum(),momenta,onShell);
// matrix element piece
wgt *= decayer->me2(-1,*inpart,outgoing_,momenta,DecayIntegrator2::Initialize);
}
catch (Veto) {
wgt=0.;
}
if(wgt>maxWeight_) maxWeight_ = wgt;
wsum += wgt;
wsqsum += sqr(wgt);
}
wsum /= decayer->nPoint_;
wsqsum=wsqsum/decayer->nPoint_-sqr(wsum);
if(wsqsum<0.) wsqsum=0.;
wsqsum=sqrt(wsqsum/decayer->nPoint_);
Energy fact = (widthGen_&&partial_>=0) ?
widthGen_->partialWidth(partial_,inpart->nominalMass()) :
inpart->dataPtr()->width();
if(fact==ZERO) fact=MeV;
// factor for the weight with spin correlations
maxWeight_ *= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6;
if ( Debug::level > 1 ) {
// ouptut the information on the initialisation
CurrentGenerator::log() << "Initialized the phase space for the decay "
<< incoming_.first->PDGName() << " -> ";
for(tPDPtr part : outgoing_)
CurrentGenerator::log() << part->PDGName() << " ";
CurrentGenerator::log() << "\n";
if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << wsum
<< " +/- " << wsqsum << "\n";
CurrentGenerator::log() << "The partial width is " << wsum*fact/MeV
<< " +/- " << wsqsum*fact/MeV << " MeV\n";
CurrentGenerator::log() << "The partial width is "
<< wsum*fact/6.58212E-22/MeV
<< " +/- " << wsqsum*fact/6.58212E-22/MeV<< " s-1\n";
CurrentGenerator::log() << "The maximum weight is "
<< maxWeight_ << endl;
}
output=wsum*fact;
}
else {
vector<Lorentz5Momentum> momenta(outgoing_.size());
double totsum(0.),totsq(0.);
for(unsigned int iy=0;iy<decayer->nIter_;++iy) {
// zero the maximum weight
maxWeight_=0.;
vector<double> wsum(channels_.size(),0.),wsqsum(channels_.size(),0.);
vector<int> nchan(channels_.size(),0);
totsum = 0.;
totsq = 0.;
Energy m0,mmin(ZERO);
for(tcPDPtr part : outgoing_) mmin += part->massMin();
for(unsigned int ix=0;ix<decayer->nPoint_;++ix) {
m0 = !onShell ? incoming_.first->generateMass() : incoming_.first->mass();
double wgt=0.;
int ichan(-1);
if(m0>mmin) {
inpart->set5Momentum(Lorentz5Momentum(m0));
// compute the prefactor
Energy prewid= (widthGen_&&partial_>=0) ?
widthGen_->partialWidth(partial_,inpart->mass()) :
inpart->dataPtr()->width();
InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
// generate the weight for this point
try {
fillStack();
wgt = pre*weight(ichan,inpart->momentum(),momenta,onShell);
// matrix element piece
wgt *= decayer->me2(-1,*inpart,outgoing_,momenta,DecayIntegrator2::Initialize);
}
catch (Veto) {
wgt=0.;
}
}
if(wgt>maxWeight_) maxWeight_=wgt;
if(ichan>=0) {
wsum[ichan] += wgt;
wsqsum[ichan] += sqr(wgt);
++nchan[ichan];
}
totsum+=wgt;
totsq+=wgt*wgt;
}
totsum=totsum/decayer->nPoint_;
totsq=totsq/decayer->nPoint_-sqr(totsum);
if(totsq<0.) totsq=0.;
totsq=sqrt(totsq/decayer->nPoint_);
if ( Debug::level > 1 )
CurrentGenerator::log() << "The branching ratio is " << iy << " "
<< totsum << " +/- " << totsq
<< maxWeight_ << "\n";
// compute the individual terms
double total(0.);
for(unsigned int ix=0;ix<channels_.size();++ix) {
if(nchan[ix]!=0) {
wsum[ix]=wsum[ix]/nchan[ix];
wsqsum[ix]=wsqsum[ix]/nchan[ix];
if(wsqsum[ix]<0.) wsqsum[ix]=0.;
wsqsum[ix]=sqrt(wsqsum[ix]/nchan[ix]);
}
else {
wsum[ix]=0;
wsqsum[ix]=0;
}
total+=sqrt(wsqsum[ix])*channels_[ix].weight();
}
if(total>0.) {
for(unsigned int ix=0;ix<channels_.size();++ix) {
channels_[ix].weight(sqrt(wsqsum[ix])*channels_[ix].weight()/total);
}
}
}
// factor for the weight with spin correlations
maxWeight_*= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6;
// output the information on the initialisation
Energy fact = (widthGen_&&partial_>=0) ?
widthGen_->partialWidth(partial_,inpart->nominalMass()) :
inpart->dataPtr()->width();
output=totsum*fact;
if(fact==ZERO) fact=MeV;
if ( Debug::level > 1 ) {
CurrentGenerator::log() << "Initialized the phase space for the decay "
<< incoming_.first->PDGName() << " -> ";
for(tcPDPtr part : outgoing_)
CurrentGenerator::log() << part->PDGName() << " ";
CurrentGenerator::log() << "\n";
if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << totsum
<< " +/- " << totsq << "\n";
CurrentGenerator::log() << "The partial width is " << totsum*fact/MeV
<< " +/- " << totsq*fact/MeV << " MeV\n";
CurrentGenerator::log() << "The partial width is "
<< totsum*fact/6.58212E-22/MeV
<< " +/- " << totsq*fact/6.58212E-22/MeV
<< " s-1\n";
CurrentGenerator::log() << "The maximum weight is "
<< maxWeight_ << "\n";
CurrentGenerator::log() << "The weights for the different phase"
<< " space channels are \n";
for(unsigned int ix=0,N=channels_.size();ix<N;++ix) {
CurrentGenerator::log() << "Channel " << ix
<< " had weight " << channels_[ix].weight()
<< "\n";
}
}
CurrentGenerator::log() << flush;
}
decayer->ME(DecayMEPtr());
return output;
}
// generate a phase-space point using multichannel phase space
Energy PhaseSpaceMode::channelPhaseSpace(int & ichan, const Lorentz5Momentum & in,
vector<Lorentz5Momentum> & momenta,
bool onShell) const {
double wgt(rStack_.top());
rStack_.pop();
// select a channel
ichan=-1;
do {
++ichan;
wgt-=channels_[ichan].weight();
}
while(ichan<int(channels_.size())&&wgt>0.);
// generate the momenta
if(ichan==int(channels_.size())) {
throw DecayIntegratorError() << "PhaseSpaceMode::channelPhaseSpace()"
<< " failed to select a channel"
<< Exception::abortnow;
}
// generate the masses of the external particles
double masswgt(1.);
vector<Energy> mass(externalMasses(in.mass(),masswgt,onShell));
momenta=channels_[ichan].generateMomenta(in,mass);
// compute the denominator of the weight
wgt=0.;
for(const PhaseSpaceChannel & channel : channels_)
wgt += channel.generateWeight(momenta);
// return the weight
return wgt!=0. ? in.mass()*masswgt/wgt : ZERO;
}
void PhaseSpaceMode::init() {
// get mass and width generators
outgoingCC_.clear();
for(tPDPtr part : outgoing_) {
outgoingCC_.push_back(part->CC() ? part->CC() : part);
}
massGen_.resize(outgoing_.size());
widthGen_ = dynamic_ptr_cast<cGenericWidthGeneratorPtr>(incoming_.first->widthGenerator());
for(unsigned int ix=0;ix<outgoing_.size();++ix) {
assert(outgoing_[ix]);
massGen_[ix]= dynamic_ptr_cast<cGenericMassGeneratorPtr>(outgoing_[ix]->massGenerator());
}
// get max energy for decays
if(!incoming_.second)
eMax_ = testOnShell_ ? incoming_.first->mass() : incoming_.first->massMax();
// work out how many random numbers we need
nRand_ = 3*outgoing_.size()-4;
for(unsigned int ix=0;ix<outgoing_.size();++ix) {
if(massGen_[ix]) ++nRand_;
}
if(channels_.empty()) return;
++nRand_;
// ensure weights sum to one
double sum(0.);
for(const PhaseSpaceChannel & channel : channels_)
sum+=channel.weight();
for(PhaseSpaceChannel & channel : channels_)
channel.weight(channel.weight()/sum);
}
void PhaseSpaceMode::initrun() {
// update the mass and width generators
if(incoming_.first->widthGenerator()!=widthGen_)
widthGen_ = dynamic_ptr_cast<cGenericWidthGeneratorPtr>(incoming_.first->widthGenerator());
for(unsigned int ix=0;ix<outgoing_.size();++ix) {
if(massGen_[ix]!=outgoing_[ix]->massGenerator())
massGen_[ix] = dynamic_ptr_cast<cGenericMassGeneratorPtr>(outgoing_[ix]->massGenerator());
}
for(PhaseSpaceChannel & channel : channels_) channel.initrun(this);
if(widthGen_) const_ptr_cast<GenericWidthGeneratorPtr>(widthGen_)->initrun();
tcGenericWidthGeneratorPtr wtemp;
for(unsigned int ix=0;ix<outgoing_.size();++ix) {
wtemp=
dynamic_ptr_cast<tcGenericWidthGeneratorPtr>(outgoing_[ix]->widthGenerator());
if(wtemp) const_ptr_cast<tGenericWidthGeneratorPtr>(wtemp)->initrun();
}
// ensure weights sum to one
double sum(0.);
for(const PhaseSpaceChannel & channel : channels_)
sum+=channel.weight();
for(PhaseSpaceChannel & channel : channels_)
channel.weight(channel.weight()/sum);
}
diff --git a/Decay/PhaseSpaceMode.h b/Decay/PhaseSpaceMode.h
--- a/Decay/PhaseSpaceMode.h
+++ b/Decay/PhaseSpaceMode.h
@@ -1,388 +1,398 @@
// -*- C++ -*-
#ifndef Herwig_PhaseSpaceMode_H
#define Herwig_PhaseSpaceMode_H
//
// This is the declaration of the PhaseSpaceMode class.
//
#include "ThePEG/Config/ThePEG.h"
#include "PhaseSpaceMode.fh"
#include "PhaseSpaceChannel.h"
#include "Herwig/PDT/GenericWidthGenerator.h"
#include "Herwig/PDT/GenericMassGenerator.h"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
*
* The <code>PhaseSpaceMode</code> class is designed to store a group
* of phase-space channels for use by the DecayIntegrator class to
* generate the phase-space for a given decay mode.
*
* Additional phase-space channels can be added using the addChannel member.
*
* In practice the modes are usually constructed together with the a number of
* <code>PhaseSpaceChannel</code> objects. In classes inheriting from the
* DecayIntegrator class.
*
* @see DecayIntegrator
* @see PhaseSpaceChannel
*
* @author Peter Richardson
*
*/
class PhaseSpaceMode: public Base {
friend class PhaseSpaceChannel;
public:
/** @name Standard constructors and destructors. */
//@{
/**
* The default constructor.
*/
PhaseSpaceMode() : maxWeight_(0.), partial_(-1),
testOnShell_(false), eMax_(ZERO),
eps_(ZERO), nRand_(0)
{};
/**
* The default constructor.
*/
PhaseSpaceMode(tPDPtr in1, tPDVector out,
double wMax, tPDPtr in2=tPDPtr(),
Energy eMax=ZERO) : incoming_(make_pair(in1,in2)),
maxWeight_(wMax),
outgoing_(out), partial_(-1),
testOnShell_(false), eMax_(eMax),
eps_(ZERO), nRand_(0)
{};
//@}
public:
/**
* Generate the decay.
* @param intermediates Whether or not to generate the intermediate particle
* in the decay channel.
* @param cc Whether we are generating the mode specified or the charge
* conjugate mode.
* @param inpart The incoming particle.
* @return The outgoing particles.
*/
ParticleVector generateDecay(const Particle & inpart,
tcDecayIntegrator2Ptr decayer,
bool intermediates,bool cc);
/**
* Add a new channel.
* @param channel A pointer to the new PhaseChannel
*/
void addChannel(PhaseSpaceChannel channel) {
channel.init(this);
channels_.push_back(channel);
}
/**
* Reset the properties of one of the intermediate particles. Only a specific channel
* is reset.
* @param ichan The channel to reset.
* @param part The ParticleData object of the particle to reset
* @param mass The mass of the intermediate.
* @param width The width of gthe intermediate.
*/
void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) {
if(!part) return;
channels_[ichan].resetIntermediate(part,mass,width);
}
/**
* Reset the properties of one of the intermediate particles. All the channels
* are reset.
* @param part The ParticleData object of the particle to reset
* @param mass The mass of the intermediate.
* @param width The width of gthe intermediate.
*/
void resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
for(PhaseSpaceChannel & channel : channels_)
channel.resetIntermediate(part,mass,width);
}
/**
* The phase-space channels
*/
const vector<PhaseSpaceChannel> & channels() const {return channels_;}
/**
* Set the weights
*/
void setWeights(const vector<double> & wgts) {
assert(wgts.size()==channels_.size());
for(unsigned int ix=0;ix<channels_.size();++ix)
channels_[ix].weight(wgts[ix]);
}
/**
+ * Access to the selected channel
+ */
+ unsigned int selectedChannel() const {return iChannel_;}
+
+ /**
* Number of rnadom numbers needed
*/
unsigned int nRand() const {return nRand_;}
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();
public :
/**
* Initialize the phase space
*/
void init();
/**
* Initialisation before the run stage
*/
void initrun();
/**
* Get the maximum weight for the decay.
* @return The maximum weight.
*/
double maxWeight() const {return maxWeight_;}
/**
* Set the maximum weight for the decay.
* @return The maximum weight.
*/
void maxWeight(double wgt) const {maxWeight_=wgt;}
/**
* Initialise the phase space.
* @param init Perform the initialization.
*/
Energy initializePhaseSpace(bool init, tcDecayIntegrator2Ptr decayer,
bool onShell=false);
/**
* The incoming particles
*/
pair<PDPtr,PDPtr> incoming() const {return incoming_;}
/**
* Access to the outging particles.
* @return A pointer to the ParticleData object.
*/
tPDVector outgoing() const {return outgoing_;}
/**
* Number of outgoing particles.
* @return The number of outgoing particles.
*/
unsigned int numberOfParticles() const {return outgoing_.size();}
/**
* Set the partial width to use for normalization. This is the partial width
* in the WidthGenerator object.
* @param in The partial width to use.
*/
void setPartialWidth(int in) {partial_=in;}
/**
* Access to the epsilon parameter
*/
Energy epsilonPS() const {return eps_;}
/**
* Fill the stack
*/
void fillStack(const double * r) {
assert(rStack_.empty());
for(unsigned int ix=nRand_;ix>0;--ix)
rStack_.push(r[nRand_-1]);
}
/**
* Fill the stack
*/
void fillStack() {
assert(rStack_.empty());
for(unsigned int ix=0;ix<nRand_;++ix) rStack_.push(UseRandom::rnd());
}
/**
* Return the weight for a given phase-space point.
* @param in The momentum of the incoming particle
* @param momenta The momenta of the outgoing particles
* @param onShell Whether or not to force the intermediates to be on-shell
* @return The weight.
*/
Energy weight(int & ichan, const Lorentz5Momentum & in,
vector<Lorentz5Momentum> & momenta,
bool onShell=false) const {
ichan=0;
// flat phase-space
if(channels_.empty())
return flatPhaseSpace(in,momenta,onShell);
// multi-channel
else
return channelPhaseSpace(ichan,in,momenta,onShell);
}
public :
/**
* A friend operator to allow the mode to be outputted for debugging purposes.
*/
friend ostream & operator<<(ostream & os, const PhaseSpaceMode & mode) {
os << "The mode has " << mode.channels_.size() << " channels\n";
if(mode.incoming_.second==PDPtr())
os << "This is a mode for the decay of " << mode.incoming_.first->PDGName() << " to ";
else
os << "This is a mode for " << mode.incoming_.first->PDGName() << ", "
<< mode.incoming_.second->PDGName() << " to ";
for(tPDPtr out : mode.outgoing_) os << out->PDGName() << " ";
os << "\n";
for(const PhaseSpaceChannel & channel : mode.channels_) os << channel;
return os;
}
private:
/**
* Return the weight and momenta for a flat phase-space decay.
* @param in The momentum of the incoming particle
* @param momenta The momenta of the outgoing particles
* @param onShell Whether or not to force the intermediates to be on-shell
* @return The weight.
*/
Energy flatPhaseSpace(const Lorentz5Momentum & in,
vector<Lorentz5Momentum> & momenta,
bool onShell=false) const;
/**
* Generate a phase-space point using multichannel phase space.
* @param ichan The channel to use
* @param in The momentum of the incoming particle
* @param momenta The momenta of the outgoing particles
* @param onShell Whether or not to force the intermediates to be on-shell
* @return The weight.
*/
Energy channelPhaseSpace(int & ichan, const Lorentz5Momentum & in,
vector<Lorentz5Momentum> & momenta,
bool onShell=false) const;
/**
* Generate the masses of the external particles.
* @param inmass The mass of the decaying particle.
* @param wgt The weight for the masses.
* @return The masses.
*/
vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;
/**
* Construct the vertex for spin corrections
* @param in The incoming particle.
* @param out The outgoing particles.
*/
void constructVertex(const Particle & in, const ParticleVector & out,
tcDecayIntegrator2Ptr decayer) const;
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
PhaseSpaceMode & operator=(const PhaseSpaceMode &);
private:
/**
* The incoming particles
*/
pair<PDPtr,PDPtr> incoming_;
/**
* The phase-space channels
*/
vector<PhaseSpaceChannel> channels_;
/**
* The maximum weight
*/
mutable double maxWeight_;
/**
* The external particles
*/
tPDVector outgoing_;
tPDVector outgoingCC_;
/**
* Which of the partial widths of the incoming particle to use
*/
int partial_;
/**
* The width generator for the incoming particle.
*/
cGenericWidthGeneratorPtr widthGen_;
/**
* The mass generators for the outgoing particles.
*/
vector<cGenericMassGeneratorPtr> massGen_;
/**
* Whether to check on-shell or off-shell kinematics
* in doinit, if on-shell off-shell is tested in initrun
*/
bool testOnShell_;
/**
* The maximum energy for the mode
*/
Energy eMax_;
/**
+ * The selected channel
+ */
+ mutable unsigned int iChannel_;
+
+ /**
* Cut-off
*/
Energy eps_;
/**
* Number of random numbers required
*/
unsigned int nRand_;
/**
* Stack for the random numbers
*/
mutable stack<double> rStack_;
};
}
#endif /* Herwig_PhaseSpaceMode_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 3:54 PM (1 d, 20 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3753976
Default Alt Text
(185 KB)

Event Timeline