Page MenuHomeHEPForge
Files F7879301

No OneTemporary
Actions

View Options

diff --git a/Decay/General/GeneralFourBodyDecayer.cc b/Decay/General/GeneralFourBodyDecayer.cc
--- a/Decay/General/GeneralFourBodyDecayer.cc
+++ b/Decay/General/GeneralFourBodyDecayer.cc
@@ -1,87 +1,918 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GeneralFourBodyDecayer class.
//
#include "GeneralFourBodyDecayer.h"
#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Switch.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"
using namespace Herwig;
-GeneralFourBodyDecayer::GeneralFourBodyDecayer() {}
-
void GeneralFourBodyDecayer::persistentOutput(PersistentOStream & os) const {
- os << _incoming << _outgoing << _diagrams << _reftag << _reftagcc;
+ os << _incoming << _outgoing << _diagrams << _reftag << _reftagcc
+ << _widthopt << _colour << _colourLargeNC << _nflow << _diagmap;
}
void GeneralFourBodyDecayer::persistentInput(PersistentIStream & is, int) {
- is >> _incoming >> _outgoing >> _diagrams >> _reftag >> _reftagcc;
+ is >> _incoming >> _outgoing >> _diagrams >> _reftag >> _reftagcc
+ >> _widthopt >> _colour >> _colourLargeNC >> _nflow >> _diagmap;
}
DescribeAbstractClass<GeneralFourBodyDecayer,DecayIntegrator>
describeGeneralFourBodyDecayer("Herwig::GeneralFourBodyDecayer",
"Herwig.so");
void GeneralFourBodyDecayer::Init() {
static ClassDocumentation<GeneralFourBodyDecayer> documentation
- ("There is no documentation for the GeneralFourBodyDecayer class");
+ ("The GeneralFourBodyDecayer class is the base class for the implementation "
+ "of all four-body decays based on spin structures in Herwig++.");
+
+ static Switch<GeneralFourBodyDecayer,unsigned int> interfaceWidthOption
+ ("WidthOption",
+ "Option for the treatment of the widths of the intermediates",
+ &GeneralFourBodyDecayer::_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);
}
ParticleVector GeneralFourBodyDecayer::decay(const Particle & parent,
const tPDVector & children) const {
- assert(false);
- return ParticleVector();
+ // 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 GeneralFourBodyDecayer::modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const {
- assert(false);
- return 0;
+ assert( !_reftag.empty() && !_reftagcc.empty() );
+ // check number of outgoing particles
+ if( children.size() != 4 || abs(parent->id()) != _incoming->id() ) return -1;
+ OrderedParticles testmode(children.begin(), children.end());
+ OrderedParticles::const_iterator dit = testmode.begin();
+ string testtag(parent->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 GeneralFourBodyDecayer::setDecayInfo(PDPtr incoming,vector<PDPtr> outgoing,
- const vector<PrototypeVertexPtr> & process) {
+
+bool GeneralFourBodyDecayer::setDecayInfo(PDPtr incoming,
+ vector<PDPtr> outgoing,
+ const vector<NBDiagram> & process) {
// set the member variables from the info supplied
// external particles
_incoming = incoming;
_outgoing = outgoing;
+ _diagrams = process;
assert( _outgoing.size() == 4 );
- // convert and store the diagrams
- for(unsigned int ix=0;ix<process.size();++ix)
- _diagrams.push_back(NBDiagram(process[ix]));
// Construct reference tags for testing in modeNumber function
OrderedParticles refmode(_outgoing.begin(), _outgoing.end());
OrderedParticles::const_iterator dit = refmode.begin();
_reftag = _incoming->name() + "->";
for( ; dit != refmode.end(); ++dit) {
if( dit != refmode.begin() ) _reftag += string(",");
_reftag += (**dit).name();
}
//CC-mode
refmode.clear();
_reftagcc = _incoming->CC() ? _incoming->CC()->name() : _incoming->name();
_reftagcc += "->";
- for( unsigned int i = 0; i < 3; ++i ) {
+ for( unsigned int i = 0; i < _outgoing.size(); ++i ) {
if( _outgoing[i]->CC() ) refmode.insert( _outgoing[i]->CC() );
else refmode.insert( _outgoing[i] );
}
dit = refmode.begin();
for( ; dit != refmode.end(); ++dit) {
if( dit != refmode.begin() ) _reftag += string(",");
_reftagcc += (**dit).name();
}
- // set the color factors
- cerr << "testing in setDecayInfo\n";
-// return setColourFactors();
- exit(0);
+ // set the colour factors and return the answer
+ return setColourFactors();
}
+
+Energy GeneralFourBodyDecayer::partialWidth(tPDPtr inpart,
+ OrderedParticles outgoing) const {
+ tPDVector temp = tPDVector(outgoing.begin(),outgoing.end());
+ bool cc=false;
+ int imode = modeNumber(cc,inpart,temp);
+ if(imode<0) return ZERO;
+ else return initializePhaseSpaceMode(0,true);
+}
+
+namespace {
+ bool addIntermediate(DecayPhaseSpaceChannelPtr newChannel,
+ const NBVertex &vertex,int &order,
+ int & ioff, int & loc) {
+ if(vertex.vertices.size()!=2) return false;
+ assert(!vertex.vertices.begin()->second.incoming);
+ int first = order/pow(10,loc);
+ order = order%int(pow(10,loc));
+ --loc;
+ // one off-shell
+ if((++vertex.vertices.begin())->second.incoming) {
+ newChannel->addIntermediate(vertex.incoming,0,0.0,ioff,first);
+ --ioff;
+ return addIntermediate(newChannel,(++vertex.vertices.begin())->second,
+ order,ioff,loc);
+ }
+ // no off-shell
+ else {
+ int second = order/pow(10,loc);
+ order = order%int(pow(10,loc));
+ --loc;
+ newChannel->addIntermediate(vertex.incoming,0,0.0,first,second);
+ --ioff;
+ return true;
+ }
+ }
+}
+
+void GeneralFourBodyDecayer::doinit() {
+ DecayIntegrator::doinit();
+ // 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;
+ // create the phase-space channels for the integration
+ unsigned int nmode(0);
+ for(unsigned int ix=0;ix<_diagrams.size();++ix) {
+ // check not four point
+ if(_diagrams[ix].vertices.size()!=2) continue;
+ // create the new channel
+ newChannel=new_ptr(DecayPhaseSpaceChannel(mode));
+ int order = _diagrams[ix].channelType;
+ int ioff = -1;
+ int loc=3;
+ // two off-shell particles
+ if(_diagrams[ix].vertices.begin()->second.incoming) {
+ newChannel->addIntermediate(extpart[0],0,0.0,-1,-2);
+ if(!addIntermediate(newChannel,( _diagrams[ix].vertices.begin())->second,
+ order,ioff,loc)) continue;
+ if(!addIntermediate(newChannel,(++_diagrams[ix].vertices.begin())->second,
+ order,ioff,loc)) continue;
+ }
+ // only one off-shell particle
+ else {
+ int first = order/pow(10,loc);
+ order = order%int(pow(10,loc));
+ --loc;
+ newChannel->addIntermediate(extpart[0],0,0.0,-1,first);
+ --ioff;
+ if(!addIntermediate(newChannel,(++_diagrams[ix].vertices.begin())->second,
+ order,ioff,loc)) continue;
+ }
+ _diagmap.push_back(ix);
+ mode->addChannel(newChannel);
+ cerr << "testing channel " << *newChannel << "\n";
+ ++nmode;
+ }
+ if(nmode==0) {
+ string mode = extpart[0]->PDGName() + "->";
+ for(unsigned int ix=1;ix<extpart.size();++ix) mode += extpart[ix]->PDGName() + " ";
+ throw Exception() << "No decay channels in GeneralFourBodyDecayer::"
+ << "doinit() for " << mode << "\n" << Exception::runerror;
+ }
+ // add the mode
+ vector<double> wgt(nmode,1./double(nmode));
+ addMode(mode,1.,wgt);
+}
+
+void GeneralFourBodyDecayer::colourConnections(const Particle & parent,
+ const ParticleVector & out) const {
+ // first extract the outgoing particles and intermediate
+ ParticleVector 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(4);
+// 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 GeneralFourBodyDecayer::"
+// << "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 GeneralFourBodyDecayer::"
+// << "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 GeneralFourBodyDecayer::"
+// << "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 "
+// << "GeneralFourBodyDecayer::"
+// << "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 "
+// << "GeneralFourBodyDecayer::"
+// << "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 {
+// string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+// + out[1]->PDGName() + " " + out[2]->PDGName();
+// throw Exception()
+// << "Unknown colour structure in GeneralFourBodyDecayer::"
+// << "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"
+// << " GeneralFourBodyDecayer::"
+// << "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 "
+// << "GeneralFourBodyDecayer::"
+// << "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 {
+// string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+// + out[1]->PDGName() + " " + out[2]->PDGName();
+// throw Exception()
+// << "Unknown colour structure in GeneralFourBodyDecayer::"
+// << "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 GeneralFourBodyDecayer::"
+// << "colourConnections() for "
+// << "decaying colour octet "
+// << mode << Exception::runerror;
+// }
+// }
+// }
+// else {
+// string mode = parent.PDGName() + " -> " + out[0]->PDGName() + " "
+// + out[1]->PDGName() + " " + out[2]->PDGName();
+// throw Exception()
+// << "Unknown colour structure in GeneralFourBodyDecayer::"
+// << "colourConnections() for octet decaying particle"
+// << mode << Exception::runerror;
+// }
+// }
+}
+
+bool GeneralFourBodyDecayer::setColourFactors() {
+ string name = _incoming->PDGName() + "->";
+ vector<int> sng,trip,atrip,oct;
+ unsigned int iloc(0);
+ for(vector<PDPtr>::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()==4) {
+ _nflow = 1;
+ _colour = vector<DVector>(1,DVector(1,1.));
+ _colourLargeNC = vector<DVector>(1,DVector(1,1.));
+ }
+ else if(sng.size()==2&&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) {
+// 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(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 getColourFactors(), omitting decay\n";
+ return false;
+ }
+ }
+ // colour triplet decaying particle
+ else if( _incoming->iColour() == PDT::Colour3) {
+ if(sng.size()==3&&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&&sng.size()==1) {
+ _nflow = 2;
+ _colour.resize(2,DVector(2,0.));
+ _colour[0][0] = 3.; _colour[0][1] = 1.;
+ _colour[1][0] = 1.; _colour[1][1] = 3.;
+ _colourLargeNC.resize(2,DVector(2,0.));
+ _colourLargeNC[0][0] = 3.; _colourLargeNC[1][1] = 3.;
+ // particle connect to incoming in first flow
+ unsigned int ifirst(0);
+ for(unsigned int ix=0;ix<_outgoing.size();++ix) {
+ if(_outgoing[ix]->iColour()==PDT::Colour3&&ifirst==0) {
+ ifirst = ix+1;
+ break;
+ }
+ }
+ // sort out the contribution of the different diagrams
+ // to the colour flows
+ for(vector<NBDiagram>::iterator it = _diagrams.begin();
+ it!=_diagrams.end();++it) {
+ // first topology
+ if(it->vertices.begin()->second.incoming) {
+ tPDPair inter = make_pair(it->vertices.begin() ->second.incoming,
+ (++(it->vertices.begin()))->second.incoming);
+ // one neutral and one triplet
+ if(inter. first->iColour()==PDT::Colour3 &&
+ inter.second->iColour()==PDT::Colour0) {
+ if( it->channelType/1000 == ifirst ||
+ (it->channelType%1000)/100 == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ else if(inter. first->iColour()==PDT::Colour0 &&
+ inter.second->iColour()==PDT::Colour3) {
+ if((it->channelType%100)/10 == ifirst ||
+ it->channelType%10 == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ // one neutral and one octet
+ else if(inter. first->iColour()==PDT::Colour3 &&
+ inter.second->iColour()==PDT::Colour8) {
+ if( it->channelType/1000 == ifirst ||
+ (it->channelType%1000)/100 == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else if(inter. first->iColour()==PDT::Colour8 &&
+ inter.second->iColour()==PDT::Colour3) {
+ if((it->channelType%100)/10 == ifirst ||
+ it->channelType%10 == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else
+ return false;
+ }
+ // second topology
+ else {
+ tPDPtr inter = (++(it->vertices.begin()))->second.incoming;
+ unsigned int iloc;
+ if(inter->iColour()==PDT::Colour3) {
+ iloc = (it->channelType%1000)/100;
+ const NBVertex & vertex = (++(it->vertices.begin()))->second;
+ inter = (++vertex.vertices.begin())->second.incoming;
+ }
+ else {
+ iloc = it->channelType/1000;
+ }
+ if(inter->iColour()==PDT::Colour0) {
+ if( iloc == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ else if(inter->iColour()==PDT::Colour8) {
+ if( iloc == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else
+ return false;
+ }
+ }
+ }
+ else {
+ generator()->log() << "Unknown colour structure for "
+ << "triplet decay in "
+ << "FourBodyDecayConstructor::getColourFactors()"
+ << " for " << name << " omitting decay\n";
+ return false;
+ }
+ }
+ // colour antitriplet decaying particle
+ else if( _incoming->iColour() == PDT::Colour3bar) {
+ if(sng.size()==3&&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&&sng.size()==1) {
+ _nflow = 2;
+ _colour.resize(2,DVector(2,0.));
+ _colour[0][0] = 3.; _colour[0][1] = 1.;
+ _colour[1][0] = 1.; _colour[1][1] = 3.;
+ _colourLargeNC.resize(2,DVector(2,0.));
+ // particle connect to incoming in first flow
+ unsigned int ifirst(0);
+ for(unsigned int ix=0;ix<_outgoing.size();++ix) {
+ if(_outgoing[ix]->iColour()==PDT::Colour3&&ifirst==0) {
+ ifirst = ix+1;
+ break;
+ }
+ }
+ // sort out the contribution of the different diagrams
+ // to the colour flows
+ for(vector<NBDiagram>::iterator it = _diagrams.begin();
+ it!=_diagrams.end();++it) {
+ // first topology
+ if(it->vertices.begin()->second.incoming) {
+ tPDPair inter = make_pair(it->vertices.begin() ->second.incoming,
+ (++(it->vertices.begin()))->second.incoming);
+ // one neutral and one triplet
+ if(inter. first->iColour()==PDT::Colour3bar &&
+ inter.second->iColour()==PDT::Colour0) {
+ if( it->channelType/1000 == ifirst ||
+ (it->channelType%1000)/100 == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ else if(inter. first->iColour()==PDT::Colour0 &&
+ inter.second->iColour()==PDT::Colour3bar) {
+ if((it->channelType%100)/10 == ifirst ||
+ it->channelType%10 == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ // one neutral and one octet
+ else if(inter. first->iColour()==PDT::Colour3bar &&
+ inter.second->iColour()==PDT::Colour8) {
+ if( it->channelType/1000 == ifirst ||
+ (it->channelType%1000)/100 == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else if(inter. first->iColour()==PDT::Colour8 &&
+ inter.second->iColour()==PDT::Colour3bar) {
+ if((it->channelType%100)/10 == ifirst ||
+ it->channelType%10 == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else
+ return false;
+ }
+ // second topology
+ else {
+ tPDPtr inter = (++(it->vertices.begin()))->second.incoming;
+ unsigned int iloc;
+ if(inter->iColour()==PDT::Colour3bar) {
+ iloc = (it->channelType%1000)/100;
+ const NBVertex & vertex = (++(it->vertices.begin()))->second;
+ inter = (++vertex.vertices.begin())->second.incoming;
+ }
+ else {
+ iloc = it->channelType/1000;
+ }
+ if(inter->iColour()==PDT::Colour0) {
+ if( iloc == ifirst) {
+ it-> colourFlow = vector<CFPair>(1,make_pair(1,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(1,1.));
+ }
+ else {
+ it->colourFlow = vector<CFPair>(1,make_pair(2,1.));
+ it->largeNcColourFlow = vector<CFPair>(1,make_pair(2,1.));
+ }
+ }
+ else if(inter->iColour()==PDT::Colour8) {
+ if( iloc == ifirst) {
+ vector<CFPair> flow(1,make_pair(2, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(1,-1./6.));
+ it->colourFlow=flow;
+ }
+ else {
+ vector<CFPair> flow(1,make_pair(1, 0.5 ));
+ it->largeNcColourFlow = flow;
+ flow.push_back( make_pair(2,-1./6.));
+ it->colourFlow=flow;
+ }
+ }
+ else
+ return false;
+ }
+ }
+ }
+ else {
+ generator()->log() << "Unknown colour antitriplet decay in "
+ << "FourBodyDecayConstructor::getColourFactors()"
+ << " for " << name << " omitting decay\n";
+ return false;
+ }
+ }
+ else if( _incoming->iColour() == PDT::Colour8) {
+ // triplet antitriplet
+ if(trip.size() == 1 && atrip.size() == 1 && sng.size() == 2) {
+ _nflow = 1;
+ _colour = vector<DVector>(1,DVector(1,0.5));
+ _colourLargeNC = vector<DVector>(1,DVector(1,0.5));
+ }
+ else {
+ generator()->log() << "Unknown colour octet decay in "
+ << "FourBodyDecayConstructor::getColourFactors()"
+ << " for " << name << " omitting decay\n";
+ return false;
+ }
+ }
+ 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";
+ }
+ }
+ return true;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Decay/General/GeneralFourBodyDecayer.h b/Decay/General/GeneralFourBodyDecayer.h
--- a/Decay/General/GeneralFourBodyDecayer.h
+++ b/Decay/General/GeneralFourBodyDecayer.h
@@ -1,159 +1,251 @@
// -*- C++ -*-
#ifndef HERWIG_GeneralFourBodyDecayer_H
#define HERWIG_GeneralFourBodyDecayer_H
//
// This is the declaration of the GeneralFourBodyDecayer class.
//
#include "Herwig++/Decay/DecayIntegrator.h"
#include "Herwig++/Models/General/PrototypeVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the GeneralFourBodyDecayer class.
*
* @see \ref GeneralFourBodyDecayerInterfaces "The interfaces"
* defined for GeneralFourBodyDecayer.
*/
class GeneralFourBodyDecayer: public DecayIntegrator {
public:
/**
* The default constructor.
*/
- GeneralFourBodyDecayer();
+ GeneralFourBodyDecayer(): _nflow(999), _widthopt(1),
+ _reftag(), _reftagcc(), _iflow(999)
+ {}
/** @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;
+ virtual int modeNumber(bool & cc, tcPDPtr parent,
+ const tPDVector & children) const;
/**
* Set the diagrams
*/
bool setDecayInfo(PDPtr incoming,vector<PDPtr> outgoing,
- const vector<PrototypeVertexPtr> & process);
+ const vector<NBDiagram> & process);
//@}
+
+ /**
+ * Function to return partial Width
+ * @param inpart Pointer to incoming particle data object
+ * @param outgoing the decay products
+ */
+ virtual Energy partialWidth(tPDPtr inpart,
+ OrderedParticles outgoing) 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:
+
+ /**
+ * Incoming particle
+ */
+ PDPtr incoming() const { return _incoming; }
+
+ /**
+ * Outgoing particles
+ */
+ const vector<PDPtr> & outgoing() const { return _outgoing; }
+
+ /**
+ * Number of colour flows
+ */
+ unsigned int numberOfFlows() const { return _nflow; }
+
+ /**
+ * Set up the colour factors
+ */
+ bool setColourFactors();
+
+ /**
+ * 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;
+ }
+
+ /**
+ * 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;
+
+ /**
+ * 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; }
+
+ /**
+ * Access the TBDiagrams that store the required information
+ * to create the diagrams
+ */
+ const vector<NBDiagram> & getProcessInfo() const {
+ return _diagrams;
+ }
+
+ /**
+ * Get the mapping between the phase-space channel and the diagram
+ */
+ const vector<unsigned int> & diagramMap() const {
+ return _diagmap;
+ }
+
+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.
*/
GeneralFourBodyDecayer & operator=(const GeneralFourBodyDecayer &);
private:
/**
* Store the incoming particle
*/
PDPtr _incoming;
/**
* Outgoing particles
*/
vector<PDPtr> _outgoing;
/**
* Store the diagrams for the decay
*/
vector<NBDiagram> _diagrams;
-// /**
-// * Map between the diagrams and the phase-space channels
-// */
-// vector<unsigned int> _diagmap;
+ /**
+ * Map between the diagrams and the phase-space channels
+ */
+ vector<unsigned int> _diagmap;
-// /**
-// * Store colour factors for ME calc.
-// */
-// vector<DVector> _colour;
+ /**
+ * Store colour factors for ME calc.
+ */
+ vector<DVector> _colour;
-// /**
-// * Store cololur factors for ME calc at large N_c
-// */
-// vector<DVector> _colourLargeNC;
+ /**
+ * Store cololur factors for ME calc at large N_c
+ */
+ vector<DVector> _colourLargeNC;
-// /**
-// * The number of colourflows.
-// */
-// unsigned int _nflow;
+ /**
+ * 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;
+ /**
+ * 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;
+ /**
+ * The colour flow
+ */
+ mutable unsigned int _iflow;
};
}
#endif /* HERWIG_GeneralFourBodyDecayer_H */
diff --git a/Decay/General/StoFFFFDecayer.cc b/Decay/General/StoFFFFDecayer.cc
--- a/Decay/General/StoFFFFDecayer.cc
+++ b/Decay/General/StoFFFFDecayer.cc
@@ -1,50 +1,409 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the StoFFFFDecayer class.
//
#include "StoFFFFDecayer.h"
#include "ThePEG/PDT/DecayMode.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++/Decay/DecayPhaseSpaceMode.h"
+#include <numeric>
using namespace Herwig;
-StoFFFFDecayer::StoFFFFDecayer() {}
-
IBPtr StoFFFFDecayer::clone() const {
return new_ptr(*this);
}
IBPtr StoFFFFDecayer::fullclone() const {
return new_ptr(*this);
}
void StoFFFFDecayer::persistentOutput(PersistentOStream & os) const {
+ os << firstVVS_ << firstVSS_ << firstSSS_ << firstFFS_
+ << secondFFV_ << secondFFS_
+ << thirdFFV_ << thirdFFS_ << sign_;
}
void StoFFFFDecayer::persistentInput(PersistentIStream & is, int) {
+ is >> firstVVS_ >> firstVSS_ >> firstSSS_ >> firstFFS_
+ >> secondFFV_ >> secondFFS_
+ >> thirdFFV_ >> thirdFFS_ >> sign_;
}
DescribeClass<StoFFFFDecayer,GeneralFourBodyDecayer>
describeStoFFFFDecayer("Herwig::StoFFFFDecayer", "Herwig.so");
void StoFFFFDecayer::Init() {
static ClassDocumentation<StoFFFFDecayer> documentation
- ("The StoFFFFDecayer class performs the 4-body decays of scalar particles in BSM models");
+ ("The StoFFFFDecayer class performs the 4-fermion "
+ "decays of scalar particles in BSM models");
}
-double StoFFFFDecayer::me2(const int ichan, const Particle & part,
+void StoFFFFDecayer::doinit() {
+ GeneralFourBodyDecayer::doinit();
+ unsigned int ndiags = getProcessInfo().size();
+ firstVVS_ .resize(ndiags);
+ firstVSS_ .resize(ndiags);
+ firstSSS_ .resize(ndiags);
+ firstFFS_ .resize(ndiags);
+ secondFFV_.resize(ndiags);
+ secondFFS_.resize(ndiags);
+ thirdFFV_ .resize(ndiags);
+ thirdFFS_ .resize(ndiags);
+ for(unsigned int ix = 0;ix < ndiags; ++ix) {
+ const NBDiagram & current = getProcessInfo()[ix];
+ // first vertex
+ firstVVS_[ix] = dynamic_ptr_cast<AbstractVVSVertexPtr>(current.vertex);
+ firstVSS_[ix] = dynamic_ptr_cast<AbstractVSSVertexPtr>(current.vertex);
+ firstSSS_[ix] = dynamic_ptr_cast<AbstractSSSVertexPtr>(current.vertex);
+ firstFFS_[ix] = dynamic_ptr_cast<AbstractFFSVertexPtr>(current.vertex);
+ // get the other vertices
+ const NBVertex & second = current.vertices.begin()->second.incoming ?
+ current.vertices.begin()->second : (++current.vertices.begin())->second;
+ // get the other vertices
+ const NBVertex & third = current.vertices.begin()->second.incoming ?
+ (++current.vertices.begin())->second : (++second.vertices.begin())->second;
+ // second vertex
+ secondFFV_[ix] = dynamic_ptr_cast<AbstractFFVVertexPtr>(second.vertex);
+ secondFFS_[ix] = dynamic_ptr_cast<AbstractFFSVertexPtr>(second.vertex);
+ // third vertex
+ thirdFFV_ [ix] = dynamic_ptr_cast<AbstractFFVVertexPtr>(third .vertex);
+ thirdFFS_ [ix] = dynamic_ptr_cast<AbstractFFSVertexPtr>(third .vertex);
+ assert( ( firstVVS_[ix] || firstVSS_[ix] ||
+ firstSSS_[ix] || firstFFS_[ix] ) &&
+ (secondFFV_[ix] || secondFFS_[ix] ) &&
+ ( thirdFFV_[ix] || thirdFFS_[ix] ));
+ // NO sign
+ switch(current.channelType) {
+ case 1234: case 1342: case 1423:
+ case 2143: case 2314: case 2431:
+ case 3124: case 3241: case 3412:
+ case 4132: case 4213: case 4321:
+ sign_.push_back( 1.);
+ break;
+ case 1243: case 1324: case 1432:
+ case 2134: case 2341: case 2413:
+ case 3142: case 3214: case 3421:
+ case 4123: case 4231: case 4312:
+ sign_.push_back(-1.);
+ break;
+ default:
+ assert(false);
+ }
+ }
+}
+
+double StoFFFFDecayer::me2(const int ichan, const Particle & inpart,
const ParticleVector & decay, MEOption meopt) const {
- assert(false);
- return 0.;
+ // particle or CC of particle
+ bool cc = (*getProcessInfo().begin()).incoming->id() != inpart.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()));
+ if(meopt==Initialize) {
+ ScalarWaveFunction::
+ calculateWaveFunctions(rho_,const_ptr_cast<tPPtr>(&inpart),
+ Helicity::incoming);
+ swave_ = ScalarWaveFunction(inpart.momentum(),inpart.dataPtr(),
+ Helicity::incoming);
+ }
+ // setup spin info when needed
+ if(meopt==Terminate) {
+ ScalarWaveFunction::
+ constructSpinInfo(const_ptr_cast<tPPtr>(&inpart),
+ Helicity::incoming,true);
+ // outgoing particles
+ for(unsigned int ix = 0; ix < 4; ++ix) {
+ SpinorWaveFunction::
+ constructSpinInfo(outwave_[ix].first,decay[ix],Helicity::outgoing,true);
+ }
+ }
+ // outgoing particles
+ for(unsigned int ix = 0; ix < 4; ++ix) {
+ SpinorWaveFunction::
+ calculateWaveFunctions(outwave_[ix].first,decay[ix],Helicity::outgoing);
+ outwave_[ix].second.resize(2);
+ if(outwave_[ix].first[0].wave().Type() == u_spinortype) {
+ 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();
+ }
+ }
+ }
+ // matrix element for the colour flows
+ vector<Complex> flows(ncf, Complex(0.)),largeflows(ncf, Complex(0.));
+ vector<DecayMatrixElement> mes(ncf,DecayMatrixElement(PDT::Spin0,
+ PDT::Spin1Half,PDT::Spin1Half,
+ PDT::Spin1Half,PDT::Spin1Half));
+ vector<DecayMatrixElement> mel(ncf,DecayMatrixElement(PDT::Spin0,
+ PDT::Spin1Half,PDT::Spin1Half,
+ PDT::Spin1Half,PDT::Spin1Half));
+ // calculate the matrix element
+ 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<NBDiagram>::const_iterator dit=getProcessInfo().begin();
+ dit!=getProcessInfo().end();++dit) {
+ if(ichan>=0&&diagramMap()[ichan]!=idiag) {
+ ++idiag;
+ continue;
+ }
+ // location of the particles
+ int iloc[4]={ dit->channelType/1000 -1,
+ (dit->channelType%1000)/100-1,
+ (dit->channelType%100)/10 -1,
+ dit->channelType%10 -1};
+ // NO sign
+ double sign = sign_[idiag];
+ // work out the type of topology
+ bool topo = dit->vertices.begin()->second.incoming;
+ const NBVertex & second = topo ?
+ dit->vertices.begin() ->second :
+ (++(dit->vertices.begin()))->second;
+ const NBVertex & third = topo ?
+ (++(dit-> vertices.begin()))->second :
+ (++(second.vertices.begin()))->second;
+ // extract the intermediates
+ tPDPair inter = make_pair(second.incoming,
+ third .incoming->CC());
+ if(cc&&inter.first ->CC()) inter.first = inter.first ->CC();
+ if(cc&&inter.second->CC()) inter.second = inter.second->CC();
+
+// //\todo remove testing
+// if(abs(inter.first ->id())!=ParticleID::t||
+// abs(inter.second->id())!=ParticleID::Wplus) {
+// ++idiag;
+// continue;
+// }
+
+ // value for the diagram
+ Complex diag(0.);
+ // first compute the last part of the diagram
+ VectorWaveFunction offVector2;
+ ScalarWaveFunction offScalar2;
+ // intermeidate scalar
+ if(inter.second->iSpin()==PDT::Spin0) {
+ if(decay[iloc[2]]->id()<0&&
+ decay[iloc[3]]->id()>0) {
+ sign *= -1.;
+ offScalar2 = thirdFFS_[idiag]->
+ evaluate(scale, widthOption(),inter.second,
+ outwave_[iloc[2]].first [ihel[iloc[2]]],
+ outwave_[iloc[3]].second[ihel[iloc[3]]]);
+ }
+ else {
+ offScalar2 = thirdFFS_[idiag]->
+ evaluate(scale, widthOption(),inter.second,
+ outwave_[iloc[3]].first [ihel[iloc[3]]],
+ outwave_[iloc[2]].second[ihel[iloc[2]]]);
+ }
+ }
+ // intermediate vector
+ else if(inter.second->iSpin()==PDT::Spin1) {
+ if(decay[iloc[2]]->id()<0&&
+ decay[iloc[3]]->id()>0) {
+ sign *= -1.;
+ offVector2 = thirdFFV_[idiag]->
+ evaluate(scale, widthOption(),inter.second,
+ outwave_[iloc[2]].first [ihel[iloc[2]]],
+ outwave_[iloc[3]].second[ihel[iloc[3]]]);
+ }
+ else {
+ offVector2 = thirdFFV_[idiag]->
+ evaluate(scale, widthOption(),inter.second,
+ outwave_[iloc[3]].first [ihel[iloc[3]]],
+ outwave_[iloc[2]].second[ihel[iloc[2]]]);
+ }
+ }
+ // unknown
+ else
+ assert(false);
+ // first topology
+ if(topo) {
+ // first intermediate
+ if(inter.first->CC()) inter.first = inter.first->CC();
+ VectorWaveFunction offVector1;
+ ScalarWaveFunction offScalar1;
+ // intermeidate scalar
+ if(inter.first->iSpin()==PDT::Spin0) {
+ if(decay[iloc[0]]->id()<0&&
+ decay[iloc[1]]->id()>0) {
+ sign *= -1.;
+ offScalar1 = secondFFS_[idiag]->
+ evaluate(scale, widthOption(),inter.first,
+ outwave_[iloc[0]].first [ihel[iloc[0]]],
+ outwave_[iloc[1]].second[ihel[iloc[1]]]);
+ }
+ else {
+ offScalar1 = secondFFS_[idiag]->
+ evaluate(scale, widthOption(),inter.first,
+ outwave_[iloc[1]].first [ihel[iloc[1]]],
+ outwave_[iloc[0]].second[ihel[iloc[0]]]);
+ }
+ }
+ // intermediate vector
+ else if(inter.first->iSpin()==PDT::Spin1) {
+ if(decay[iloc[0]]->id()<0&&
+ decay[iloc[1]]->id()>0) {
+ sign *= -1.;
+ offVector1 = secondFFV_[idiag]->
+ evaluate(scale, widthOption(),inter.first,
+ outwave_[iloc[0]].first [ihel[iloc[0]]],
+ outwave_[iloc[1]].second[ihel[iloc[1]]]);
+ }
+ else {
+ offVector1 = secondFFV_[idiag]->
+ evaluate(scale, widthOption(),inter.first,
+ outwave_[iloc[1]].first [ihel[iloc[1]]],
+ outwave_[iloc[0]].second[ihel[iloc[0]]]);
+ }
+ }
+ // unknown
+ else
+ assert(false);
+ // put it all together
+ if(inter.first->iSpin()==PDT::Spin0) {
+ if(inter.second->iSpin()==PDT::Spin0) {
+ diag = firstSSS_[idiag]->
+ evaluate(scale,swave_,offScalar1,offScalar2);
+ }
+ else if(inter.second->iSpin()==PDT::Spin1) {
+ diag = firstVSS_[idiag]->
+ evaluate(scale,offVector2,offScalar1,swave_);
+ }
+ }
+ else if(inter.first->iSpin()==PDT::Spin1) {
+ if(inter.second->iSpin()==PDT::Spin0) {
+ diag = firstVSS_[idiag]->
+ evaluate(scale,offVector1,offScalar2,swave_);
+ }
+ else if(inter.second->iSpin()==PDT::Spin1) {
+ diag = firstVVS_[idiag]->
+ evaluate(scale,offVector1,offVector2,swave_);
+ }
+ }
+ }
+ // second topology
+ else {
+ if(decay[iloc[0]]->id()<0&&decay[iloc[1]]->id()>0) {
+ sign *= -1.;
+ SpinorWaveFunction inters = firstFFS_[idiag]->
+ evaluate(scale,widthOption(),inter.first,
+ outwave_[iloc[0]].first [ihel[iloc[0]]],swave_);
+ if(inter.second->iSpin()==PDT::Spin0) {
+ diag = secondFFS_[idiag]->
+ evaluate(scale,inters,outwave_[iloc[1]].second[ihel[iloc[1]]],
+ offScalar2);
+ }
+ else {
+ diag = secondFFV_[idiag]->
+ evaluate(scale,inters,outwave_[iloc[1]].second[ihel[iloc[1]]],
+ offVector2);
+ }
+ }
+ else {
+ SpinorBarWaveFunction inters = firstFFS_[idiag]->
+ evaluate(scale,widthOption(),inter.first,
+ outwave_[iloc[0]].second[ihel[iloc[0]]],swave_);
+ if(inter.second->iSpin()==PDT::Spin0) {
+ diag = secondFFS_[idiag]->
+ evaluate(scale,outwave_[iloc[1]].first [ihel[iloc[1]]],inters,
+ offScalar2);
+ }
+ else {
+ diag = secondFFV_[idiag]->
+ evaluate(scale,outwave_[iloc[1]].first [ihel[iloc[1]]],inters,
+ offVector2);
+ }
+ }
+ }
+ // 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](0,ihel[0],ihel[1],ihel[2],ihel[3]) = flows[ix];
+ mel[ix](0,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*scale*UnitRemoval::InvE2;
}
diff --git a/Decay/General/StoFFFFDecayer.h b/Decay/General/StoFFFFDecayer.h
--- a/Decay/General/StoFFFFDecayer.h
+++ b/Decay/General/StoFFFFDecayer.h
@@ -1,100 +1,184 @@
// -*- C++ -*-
#ifndef HERWIG_StoFFFFDecayer_H
#define HERWIG_StoFFFFDecayer_H
//
// This is the declaration of the StoFFFFDecayer class.
//
#include "GeneralFourBodyDecayer.h"
+#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
+#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.h"
namespace Herwig {
using namespace ThePEG;
/**
* Here is the documentation of the StoFFFFDecayer class.
*
* @see \ref StoFFFFDecayerInterfaces "The interfaces"
* defined for StoFFFFDecayer.
*/
class StoFFFFDecayer: public GeneralFourBodyDecayer {
public:
/**
- * The default constructor.
- */
- StoFFFFDecayer();
-
- /**
* 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 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;
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:
-// If needed, insert declarations of virtual function defined in the
-// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
-
+ /** @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.
*/
StoFFFFDecayer & operator=(const StoFFFFDecayer &);
+private:
+
+ /**
+ * Signs for NO
+ */
+ vector<double> sign_;
+
+ /**
+ * Potential vertices for the first step
+ */
+ //@{
+ /**
+ * VVS Vertex
+ */
+ vector<AbstractVVSVertexPtr> firstVVS_;
+
+ /**
+ * VSS Vertex
+ */
+ vector<AbstractVSSVertexPtr> firstVSS_;
+
+ /**
+ * SSS Vertex
+ */
+ vector<AbstractSSSVertexPtr> firstSSS_;
+
+ /**
+ * FFS Vertex
+ */
+ vector<AbstractFFSVertexPtr> firstFFS_;
+ //@}
+
+ /**
+ * Potential vertices for the second step
+ */
+ //@{
+ /**
+ * FFV Vertex
+ */
+ vector<AbstractFFVVertexPtr> secondFFV_;
+
+ /**
+ * FFS Vertex
+ */
+ vector<AbstractFFSVertexPtr> secondFFS_;
+ //@}
+
+ /**
+ * Potential vertices for the third step
+ */
+ //@{
+ /**
+ * FFV Vertex
+ */
+ vector<AbstractFFVVertexPtr> thirdFFV_;
+
+ /**
+ * FFS Vertex
+ */
+ vector<AbstractFFSVertexPtr> thirdFFS_;
+ //@}
+
+ /**
+ * Spin density matrix
+ */
+ mutable RhoDMatrix rho_;
+
+ /**
+ * Scalar wavefunction
+ */
+ mutable ScalarWaveFunction swave_;
+
+ /**
+ * Spinors for outgoing particles
+ */
+ mutable pair<vector<SpinorWaveFunction>,vector<SpinorBarWaveFunction> > outwave_[4];
+
};
}
#endif /* HERWIG_StoFFFFDecayer_H */
diff --git a/Models/General/FourBodyDecayConstructor.cc b/Models/General/FourBodyDecayConstructor.cc
--- a/Models/General/FourBodyDecayConstructor.cc
+++ b/Models/General/FourBodyDecayConstructor.cc
@@ -1,232 +1,258 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FourBodyDecayConstructor class.
//
#include "FourBodyDecayConstructor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/RefVector.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/PDT/DecayMode.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "Herwig++/Decay/General/GeneralFourBodyDecayer.h"
#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
#include "DecayConstructor.h"
#include <queue>
using namespace Herwig;
FourBodyDecayConstructor::~FourBodyDecayConstructor() {}
IBPtr FourBodyDecayConstructor::clone() const {
return new_ptr(*this);
}
IBPtr FourBodyDecayConstructor::fullclone() const {
return new_ptr(*this);
}
void FourBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
- os << interopt_ << widthopt_;
+ os << interOpt_ << widthOpt_;
}
void FourBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
- is >> interopt_ >> widthopt_;
+ is >> interOpt_ >> widthOpt_;
}
DescribeClass<FourBodyDecayConstructor,NBodyDecayConstructorBase>
describeFourBodyDecayConstructor("Herwig::FourBodyDecayConstructor","Herwig.so");
void FourBodyDecayConstructor::Init() {
static ClassDocumentation<FourBodyDecayConstructor> documentation
("The FourBodyDecayConstructor class implements a small number"
" of 4-body decays in general models");
static Switch<FourBodyDecayConstructor,unsigned int> interfaceWidthOption
("WidthOption",
"Option for the treatment of the widths of the intermediates",
- &FourBodyDecayConstructor::widthopt_, 1, false, false);
+ &FourBodyDecayConstructor::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<FourBodyDecayConstructor,unsigned int> interfaceIntermediateOption
("IntermediateOption",
"Option for the inclusion of intermediates in the event",
- &FourBodyDecayConstructor::interopt_, 0, false, false);
+ &FourBodyDecayConstructor::interOpt_, 0, false, false);
static SwitchOption interfaceIntermediateOptionAlways
(interfaceIntermediateOption,
"Always",
"Always include the intermediates",
1);
static SwitchOption interfaceIntermediateOptionNever
(interfaceIntermediateOption,
"Never",
"Never include the intermediates",
2);
static SwitchOption interfaceIntermediateOptionOnlyIfOnShell
(interfaceIntermediateOption,
"OnlyIfOnShell",
"Only if there are on-shell diagrams",
0);
}
void FourBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
if( particles.empty() ) return;
NBodyDecayConstructorBase::DecayList(particles);
}
void FourBodyDecayConstructor::
createDecayMode(vector<PrototypeVertexPtr> & diagrams) {
- // incoming particle
- tPDPtr inpart = diagrams[0]->incoming;
// some basic checks for the modes we are interested in
// only looking at scalars
if(diagrams[0]->incoming->iSpin()!=PDT::Spin0) return;
// which decay to 4 fermions
unsigned int nferm=0;
for(OrderedParticles::const_iterator it=diagrams[0]->outPart.begin();
it!=diagrams[0]->outPart.end();++it) {
if((**it).iSpin()==PDT::Spin1Half) ++nferm;
}
if(nferm!=4) return;
- cerr << "!!!!!!!!!!!!!!!!!! MODE !!!!!!!!!!!\n";
- cerr << "Number of diagrams " << diagrams.size() << "\n";
- cerr << diagrams[0]->incoming->PDGName() << " -> ";
- for(OrderedParticles::const_iterator it=diagrams[0]->outPart.begin();
- it!=diagrams[0]->outPart.end();++it)
- cerr << (**it).PDGName() << " ";
- cerr << "\n";
+ // check for on-shell intermediates
+ bool possibleOnShell=false;
for(unsigned int iy=0;iy<diagrams.size();++iy)
- cerr << "DIAGRAM " << iy << "\n" << *diagrams[iy] << "\n";
-// // outgoing particles
-// OrderedParticles outgoing=diagrams[0]->outPart;
-// // incoming particle is now unstable
-// inpart->stable(false);
-// // construct the tag for the decay mode
-// string tag = inpart->name() + "->";
-// for(OrderedParticles::const_iterator it = outgoing.begin();
-// it != outgoing.end(); ++it) {
-// if(it!=outgoing.begin()) tag += ",";
-// tag += (**it).name();
-// }
-// tag += ";";
-// tDMPtr dm = generator()->findDecayMode(tag);
-// if( decayConstructor()->disableDecayMode(tag) ) {
-// // If mode alread exists, ie has been read from file,
-// // disable it
-// if( dm ) {
-// generator()->preinitInterface(dm, "BranchingRatio", "set", "0.0");
-// generator()->preinitInterface(dm, "OnOff", "set", "Off");
-// }
-// return;
-// }
-// if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
-// cerr << "testing calling create " << tag << "\n";
-// // create the decayer
-// GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter);
-// if(!decayer) return;
-// }
-// else if (dm) {
-// assert(false);
-// }
-// //update CC mode if it exists
-// if( inpart->CC() )
-// inpart->CC()->synchronize();
+ possibleOnShell |= diagrams[iy]->possibleOnShell;
+ bool inter = interOpt_ == 1 || (interOpt_ == 0 && possibleOnShell);
+ // incoming particle
+ tPDPtr inpart = diagrams[0]->incoming;
+ // outgoing particles
+ OrderedParticles outgoing=diagrams[0]->outPart;
+ // incoming particle is now unstable
+ inpart->stable(false);
+ // construct the tag for the decay mode
+ string tag = inpart->name() + "->";
+ for(OrderedParticles::const_iterator it = outgoing.begin();
+ it != outgoing.end(); ++it) {
+ if(it!=outgoing.begin()) tag += ",";
+ tag += (**it).name();
+ }
+ tag += ";";
+ tDMPtr dm = generator()->findDecayMode(tag);
+ // if mode disabled zero BR and return
+ if( decayConstructor()->disableDecayMode(tag) ) {
+ // If mode alread exists, ie has been read from file,
+ // disable it
+ if( dm ) {
+ generator()->preinitInterface(dm, "BranchingRatio", "set", "0.0");
+ generator()->preinitInterface(dm, "OnOff", "set", "Off");
+ }
+ return;
+ }
+ // create mode if needed
+ if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
+ // create the decayer
+ GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter);
+ if(!decayer) {
+ if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
+ << tag << " so mode not created\n";
+ return;
+ }
+ // create the decay mode
+ tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
+ if(ndm) {
+ generator()->preinitInterface(ndm, "Decayer", "set",
+ decayer->fullName());
+ generator()->preinitInterface(ndm, "OnOff", "set", "On");
+ Energy width =
+ decayer->partialWidth(inpart,outgoing);
+ setBranchingRatio(ndm, width);
+ }
+ else
+ throw NBodyDecayConstructorError()
+ << "FourBodyDecayConstructor::createDecayMode - Needed to create "
+ << "new decaymode but one could not be created for the tag "
+ << tag << Exception::warning;
+ }
+ // otherwise
+ else if (dm && (dm->decayer()->fullName()).find("Mambo") != string::npos) {
+ // create the decayer
+ GeneralFourBodyDecayerPtr decayer = createDecayer(diagrams,inter);
+ if(!decayer) {
+ if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
+ << tag << " so mode not created\n";
+ return;
+ }
+ generator()->preinitInterface(dm, "Decayer", "set",
+ decayer->fullName());
+ }
+ //update CC mode if it exists
+ if( inpart->CC() )
+ inpart->CC()->synchronize();
}
GeneralFourBodyDecayerPtr
FourBodyDecayConstructor::createDecayer(vector<PrototypeVertexPtr> & diagrams,
bool inter) const {
if(diagrams.empty()) return GeneralFourBodyDecayerPtr();
// extract the external particles for the process
PDPtr incoming = diagrams[0]->incoming;
// outgoing particles
vector<PDPtr> outgoing(diagrams[0]->outPart.begin(),
diagrams[0]->outPart.end());
+ vector<NBDiagram> newDiagrams;
+ cerr << "!!!!!!!!!!!!!!!!!! MODE !!!!!!!!!!!\n";
+ cerr << "Number of diagrams " << diagrams.size() << "\n";
+ cerr << diagrams[0]->incoming->PDGName() << " -> ";
+ for(OrderedParticles::const_iterator it=diagrams[0]->outPart.begin();
+ it!=diagrams[0]->outPart.end();++it)
+ cerr << (**it).PDGName() << " ";
+ cerr << "\n";
+ // convert the diagrams
+ for(unsigned int ix=0;ix<diagrams.size();++ix) {
+ cerr << "DIAGRAM " << ix << "\n" << *diagrams[ix] << "\n";
+ newDiagrams.push_back(NBDiagram(diagrams[ix]));
+ }
// get the name for the object
string objectname ("/Herwig/Decays/");
string classname = DecayerClassName(incoming, diagrams[0]->outPart, objectname);
if(classname=="") return GeneralFourBodyDecayerPtr();
// create the object
GeneralFourBodyDecayerPtr decayer =
dynamic_ptr_cast<GeneralFourBodyDecayerPtr>
(generator()->preinitCreate(classname, objectname));
- // set up the decayer
- cerr << "testing made the decayer\n";
- decayer->setDecayInfo(incoming,outgoing,diagrams);
-
-
-// // get the colour flows
-// unsigned int ncf(0);
-// pair<vector<DVector>, vector<DVector> > cfactors;
-// try {
-// cfactors = getColourFactors(incoming,outgoing,diagrams,ncf);
-// }
-// catch ( Veto ) { return GeneralThreeBodyDecayerPtr(); }
-// // set decayer options from base class
-// setDecayerInterfaces(objectname);
-// // set the width option
-// ostringstream value;
-// value << _widthopt;
-// generator()->preinitInterface(objectname, "WidthOption", "set", value.str());
-// // set the intermediates option
-// ostringstream value2;
-// value2 << inter;
-// generator()->preinitInterface(objectname, "GenerateIntermediates", "set",
-// value2.str());
-// // initialize the decayer
-// decayer->init();
-// // return the decayer
-// return decayer;
-
-
- assert(false);
+ // set up the decayer and return if doesn't work
+ if(!decayer->setDecayInfo(incoming,outgoing,newDiagrams))
+ return GeneralFourBodyDecayerPtr();
+ // set decayer options from base class
+ setDecayerInterfaces(objectname);
+ // set the width option
+ ostringstream value;
+ value << widthOpt_;
+ generator()->preinitInterface(objectname, "WidthOption", "set", value.str());
+ // set the intermediates option
+ ostringstream value2;
+ value2 << inter;
+ generator()->preinitInterface(objectname, "GenerateIntermediates", "set",
+ value2.str());
+ // initialize the decayer
+ decayer->init();
+ // return the decayer
+ return decayer;
}
string FourBodyDecayConstructor::DecayerClassName(tcPDPtr incoming,
const OrderedParticles & outgoing,
string & objname) const {
string classname("Herwig::");
// spins of the outgoing particles
unsigned int ns(0),nf(0),nv(0);
objname += incoming->PDGName() + "2";
for(OrderedParticles::const_iterator it=outgoing.begin();
it!=outgoing.end();++it) {
if ((**it).iSpin()==PDT::Spin0 ) ++ns;
else if((**it).iSpin()==PDT::Spin1Half) ++nf;
else if((**it).iSpin()==PDT::Spin1 ) ++nv;
objname += (**it).PDGName();
}
objname += "Decayer";
if(incoming->iSpin()==PDT::Spin0) {
if(nf==4) classname += "StoFFFFDecayer";
else classname = "";
}
else {
classname="";
}
return classname;
}
diff --git a/Models/General/FourBodyDecayConstructor.h b/Models/General/FourBodyDecayConstructor.h
--- a/Models/General/FourBodyDecayConstructor.h
+++ b/Models/General/FourBodyDecayConstructor.h
@@ -1,141 +1,141 @@
// -*- C++ -*-
#ifndef THEPEG_FourBodyDecayConstructor_H
#define THEPEG_FourBodyDecayConstructor_H
//
// This is the declaration of the FourBodyDecayConstructor class.
//
#include "NBodyDecayConstructorBase.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "Herwig++/Decay/General/GeneralFourBodyDecayer.fh"
#include "PrototypeVertex.h"
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
/**
* Here is the documentation of the FourBodyDecayConstructor class.
*
* @see \ref FourBodyDecayConstructorInterfaces "The interfaces"
* defined for FourBodyDecayConstructor.
*/
class FourBodyDecayConstructor: public NBodyDecayConstructorBase {
public:
/**
* The default constructor.
*/
FourBodyDecayConstructor() :
- interopt_(0), widthopt_(1) {}
+ interOpt_(0), widthOpt_(1) {}
/**
* Destructor
*/
~FourBodyDecayConstructor();
/**
* Function used to determine allowed decaymodes, to be implemented
* in derived class.
* @param particles vector of ParticleData pointers containing
* particles in model
*/
virtual void DecayList(const set<PDPtr> & particles);
/**
* Number of outgoing lines. Required for correct ordering.
*/
virtual unsigned int numBodies() const {return 4;}
/**
* Create a decay mode
*/
void createDecayMode(vector<PrototypeVertexPtr> &);
/**
* Create the decayer
* @param diagrams The diagrams for the decay
* @param inter Option for intermediates
*/
GeneralFourBodyDecayerPtr createDecayer(vector<PrototypeVertexPtr> & diagrams,
bool inter) const;
/**
* Contruct the classname and object name for the Decayer
* @param incoming The incoming particle
* @param outgoing The decay products
* @param objname a string containing the default path of the Decayer object
*/
string DecayerClassName(tcPDPtr incoming, const OrderedParticles & outgoing,
string & objname) 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;
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
FourBodyDecayConstructor & operator=(const FourBodyDecayConstructor &);
private:
/**
* Option for the inclusion of intermediates
*/
- unsigned int interopt_;
+ unsigned int interOpt_;
/**
* How to treat the widths of the intermediate particles
*/
- unsigned int widthopt_;
+ unsigned int widthOpt_;
};
}
#endif /* THEPEG_FourBodyDecayConstructor_H */
diff --git a/Models/General/PrototypeVertex.cc b/Models/General/PrototypeVertex.cc
--- a/Models/General/PrototypeVertex.cc
+++ b/Models/General/PrototypeVertex.cc
@@ -1,251 +1,297 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the FourBodyDecayConstructor class.
//
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Utilities/Exception.h"
#include "PrototypeVertex.h"
using namespace Herwig;
void PrototypeVertex::createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
std::queue<PrototypeVertexPtr> & prototypes) {
int id = inpart->id();
if(!vertex->isIncoming(inpart)) return;
for(unsigned int list=0;list<vertex->getNpoint();++list) {
tPDVector decaylist = vertex->search(list, inpart);
tPDVector::size_type nd = decaylist.size();
for( tPDVector::size_type i = 0; i < nd; i += vertex->getNpoint() ) {
tPDVector pout(decaylist.begin()+i,
decaylist.begin()+i+vertex->getNpoint());
OrderedVertices out;
for(unsigned int ix=1;ix<pout.size();++ix) {
if(pout[ix]->id() == id ) swap(pout[0], pout[ix]);
if(pout[ix]->CC()) pout[ix] = pout[ix]->CC();
out.insert(make_pair(pout[ix],PrototypeVertexPtr()));
}
if(vertex->getNpoint()==3) {
// remove radiation
if((pout[0]==pout[1] && (pout[2]->id()==ParticleID::gamma||
pout[2]->id()==ParticleID::g||
pout[2]->id()==ParticleID::Z0)) ||
(pout[0]==pout[2] && (pout[1]->id()==ParticleID::gamma||
pout[1]->id()==ParticleID::g||
pout[1]->id()==ParticleID::Z0)))
continue;
}
prototypes.push(new_ptr(PrototypeVertex(inpart,out,vertex,
int(vertex->getNpoint())-1)));
}
}
}
PrototypeVertexPtr PrototypeVertex::replicateTree(PrototypeVertexPtr parent,
PrototypeVertexPtr oldChild,
PrototypeVertexPtr & newChild) {
PrototypeVertexPtr newParent =
new_ptr(PrototypeVertex(parent->incoming,OrderedVertices(),
parent->vertex,parent->npart));
for(OrderedVertices::const_iterator it = parent->outgoing.begin();
it!=parent->outgoing.end();++it) {
PrototypeVertexPtr child = it->second ?
replicateTree(it->second,oldChild,newChild) :
PrototypeVertexPtr();
newParent->outgoing.insert(make_pair(it->first,child));
if(child) child->parent = newParent;
if(it->second==oldChild) newChild=child;
}
if(oldChild==parent) newChild=newParent;
return newParent;
}
void PrototypeVertex::expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
std::queue<PrototypeVertexPtr> & prototypes,
const set<PDPtr> & excluded) {
for(OrderedVertices::const_iterator it = proto->outgoing.begin();
it!=proto->outgoing.end();++it) {
if(it->second) {
expandPrototypes(it->second,vertex,prototypes,excluded);
}
else {
if(!vertex->isIncoming(it->first)) continue;
if(excluded.find(it->first)!=excluded.end()) continue;
if(it->first->CC() &&
excluded.find(it->first->CC())!=excluded.end()) continue;
int id = it->first->id();
PrototypeVertexPtr parent=proto;
while(parent->parent) parent=parent->parent;
for(unsigned int il = 0; il < vertex->getNpoint(); ++il) {
tPDVector decaylist = vertex->search(il,it->first );
tPDVector::size_type nd = decaylist.size();
for( tPDVector::size_type i = 0; i < nd; i += vertex->getNpoint() ) {
tPDVector pout(decaylist.begin()+i,
decaylist.begin()+i+vertex->getNpoint());
OrderedVertices outgoing;
for(unsigned int iy=1;iy<pout.size();++iy) {
if(pout[iy]->id() == id ) swap(pout[0], pout[iy]);
if(pout[iy]->CC()) pout[iy] = pout[iy]->CC();
outgoing.insert(make_pair(pout[iy],PrototypeVertexPtr()));
}
if(vertex->getNpoint()==3) {
if((pout[0]==pout[1] && (pout[2]->id()==ParticleID::gamma||
pout[2]->id()==ParticleID::g||
pout[2]->id()==ParticleID::Z0)) ||
(pout[0]==pout[2] && (pout[1]->id()==ParticleID::gamma||
pout[1]->id()==ParticleID::g||
pout[1]->id()==ParticleID::Z0)))
continue;
// remove weak decays of quarks other than top
if(StandardQCDPartonMatcher::Check(pout[0]->id()) &&
((StandardQCDPartonMatcher::Check(pout[1]->id()) &&
abs(pout[2]->id())==ParticleID::Wplus)||
(StandardQCDPartonMatcher::Check(pout[2]->id())&&
abs(pout[1]->id())==ParticleID::Wplus))) continue;
// remove weak decays of leptons
if((abs(pout[0]->id())>=11&&abs(pout[0]->id())<=16) &&
(((abs(pout[1]->id())>=11&&abs(pout[1]->id())<=16) &&
abs(pout[2]->id())==ParticleID::Wplus)||
((abs(pout[2]->id())>=11&&abs(pout[2]->id())<=16)&&
abs(pout[1]->id())==ParticleID::Wplus))) continue;
}
PrototypeVertexPtr newBranch =
new_ptr(PrototypeVertex(it->first,
outgoing,vertex,int(vertex->getNpoint())-1));
PrototypeVertexPtr newChild;
PrototypeVertexPtr newVertex = replicateTree(parent,proto,newChild);
newBranch->parent = newChild;
OrderedVertices::iterator kt = newChild->outgoing.begin();
for(OrderedVertices::const_iterator jt = proto->outgoing.begin();
jt!=it;++jt,++kt) {;}
pair< tPDPtr, PrototypeVertexPtr > newPair = make_pair(kt->first,newBranch);
newChild->outgoing.erase(kt);
newChild->outgoing.insert(newPair);
newChild->incrementN(newBranch->npart-1);
prototypes.push(newVertex);
}
}
}
}
}
bool PrototypeVertex::canBeOnShell(unsigned int opt,Energy maxMass,bool first) {
Energy outMass = outgoingMass();
if(!first) {
bool in = maxMass>incomingMass();
bool out = incomingMass()>outMass;
if(opt!=0) {
if( in && ( out || opt==2 ) ) return true;
}
else if (incoming->width() == ZERO) {
tPrototypeVertexPtr original = this;
while(original->parent) {
original=original->parent;
};
ostringstream name;
name << original->incoming->PDGName() << " -> ";
for(OrderedParticles::const_iterator it = original->outPart.begin();
it!= original->outPart.end();++it)
name << (**it).PDGName() << " ";
Throw<InitException>()
<< "Trying to include on-shell diagram in decay"
<< name.str() << "including on-shell "
<< incoming->PDGName() << " with zero width.\n"
<< "You should make sure that the width for the intermediate is either"
<< " read from an SLHA file or the intermediate is included in the "
<< "DecayParticles list of the ModelGenerator.\n"
<< Exception::runerror;
}
}
else maxMass = incomingMass();
// check the decay products
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
if(!it->second) continue;
Energy newMax = maxMass-outMass+it->second->outgoingMass();
if(it->second->canBeOnShell(opt,newMax,false)) {
if(first) possibleOnShell = true;
return true;
}
}
return false;
}
+namespace {
+ void constructTag(unsigned int & loc,unsigned int & order,NBVertex & vertex,
+ const OrderedParticles & outgoing,
+ vector<bool> & matched) {
+ for(list<pair<PDPtr,NBVertex> >::iterator it=vertex.vertices.begin();
+ it!=vertex.vertices.end();++it) {
+ // search down the tree
+ if(it->second.vertex) {
+ constructTag(loc,order,it->second,outgoing,matched);
+ }
+ // identify this particle
+ else {
+ unsigned int ix=0;
+ for(OrderedParticles::const_iterator pit=outgoing.begin();
+ pit!=outgoing.end();++pit) {
+ if(*pit==it->first&&!matched[ix]) {
+ matched[ix] = true;
+ order += (ix+1)*pow(10,loc);
+ --loc;
+ break;
+ }
+ ++ix;
+ }
+ }
+ }
+ }
+}
+
NBDiagram::NBDiagram(PrototypeVertexPtr proto)
- : incoming(proto->incoming), outgoing(proto->outPart),
- vertex(proto->vertex) {
+ : channelType(0),
+ colourFlow (vector<CFPair>(1,make_pair(1,1.))),
+ largeNcColourFlow(vector<CFPair>(1,make_pair(1,1.))) {
+ if(!proto) return;
+ incoming = proto->incoming;
+ outgoing = proto->outPart;
+ vertex = proto->vertex;
// create the vertices
for(OrderedVertices::const_iterator it=proto->outgoing.begin();
it!=proto->outgoing.end();++it) {
vertices.push_back(make_pair(it->first,NBVertex(it->second)));
}
// now let's re-order so that branchings are at the end
for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
it!=vertices.end();++it) {
if(!it->second.incoming) continue;
list<pair<PDPtr,NBVertex> >::iterator jt=it;
for( ; jt!=vertices.end();++jt) {
if(jt==it) continue;
if(!jt->second.incoming) {
break;
}
}
if(jt!=vertices.end()) {
list<pair<PDPtr,NBVertex> >::iterator kt = it;
while(kt!=jt) {
list<pair<PDPtr,NBVertex> >::iterator lt = kt;
++lt;
swap(*kt,*lt);
kt=lt;
}
}
}
// finally work out the channel and the ordering
- unsigned int loc(outgoing.size());
+ unsigned int order(0);
+ unsigned int loc(outgoing.size()-1);
+ vector<bool> matched(outgoing.size(),false);
for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
it!=vertices.end();++it) {
- unsigned int ipart(0);
- for( ; ipart<outgoing.size(); ++ipart) {
+ // search down the tree
+ if(it->second.vertex) {
+ constructTag(loc,order,it->second,outgoing,matched);
}
-
-
- --loc;
+ // identify this particle
+ else {
+ unsigned int ix=0;
+ for(OrderedParticles::const_iterator pit=outgoing.begin();
+ pit!=outgoing.end();++pit) {
+ if(*pit==it->first&&!matched[ix]) {
+ matched[ix] = true;
+ order += (ix+1)*pow(10,loc);
+ --loc;
+ break;
+ }
+ ++ix;
+ }
+ }
}
-
- exit(0);
+ channelType=order;
}
NBVertex::NBVertex(PrototypeVertexPtr proto) {
if(!proto) return;
incoming = proto->incoming;
outgoing = proto->outPart;
vertex = proto->vertex;
for(OrderedVertices::const_iterator it=proto->outgoing.begin();
it!=proto->outgoing.end();++it) {
vertices.push_back(make_pair(it->first,NBVertex(it->second)));
}
// now let's re-order so that branchings are at the end
for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
it!=vertices.end();++it) {
if(!it->second.incoming) continue;
list<pair<PDPtr,NBVertex> >::iterator jt=it;
for( ; jt!=vertices.end();++jt) {
if(jt==it) continue;
if(!jt->second.incoming) {
break;
}
}
if(jt!=vertices.end()) {
list<pair<PDPtr,NBVertex> >::iterator kt = it;
while(kt!=jt) {
list<pair<PDPtr,NBVertex> >::iterator lt = kt;
++lt;
swap(*kt,*lt);
kt=lt;
}
}
}
}
diff --git a/Models/General/PrototypeVertex.h b/Models/General/PrototypeVertex.h
--- a/Models/General/PrototypeVertex.h
+++ b/Models/General/PrototypeVertex.h
@@ -1,440 +1,455 @@
// -*- C++ -*-
//
// PrototypeVertex.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2007 The Herwig Collaboration
//
// Herwig++ is licenced under version 2 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_PrototypeVertex_H
#define HERWIG_PrototypeVertex_H
#include <queue>
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/EnumIO.h"
//
// This is the declaration of the PrototypeVertex class.
//
namespace Herwig {
using namespace ThePEG;
using Helicity::VertexBasePtr;
class PrototypeVertex;
ThePEG_DECLARE_POINTERS(Herwig::PrototypeVertex,PrototypeVertexPtr);
/** Pair of int,double */
typedef pair<unsigned int, double> CFPair;
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct ParticleOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
- bool operator()(PDPtr p1, PDPtr p2) {
+ bool operator() (PDPtr p1, PDPtr p2) {
return abs(p1->id()) > abs(p2->id()) ||
( abs(p1->id()) == abs(p2->id()) && p1->id() > p2->id() ) ||
( p1->id() == p2->id() && p1->fullName() > p2->fullName() );
}
};
/**
* A struct to order the particles in the same way as in the DecayMode's
*/
struct VertexOrdering {
/**
* Operator for the ordering
* @param p1 The first ParticleData object
* @param p2 The second ParticleData object
*/
bool operator()(pair< tPDPtr, PrototypeVertexPtr > p1,
pair< tPDPtr, PrototypeVertexPtr > p2) {
return abs(p1.first->id()) > abs(p2.first->id()) ||
( abs(p1.first->id()) == abs(p2.first->id()) && p1.first->id() > p2.first->id() ) ||
( p1.first->id() == p2.first->id() && p1.first->fullName() > p2.first->fullName() );
}
};
typedef multiset<pair< tPDPtr, PrototypeVertexPtr >,VertexOrdering > OrderedVertices;
/**
* A set of ParticleData objects ordered as for the DecayMode's
*/
typedef multiset<PDPtr,ParticleOrdering> OrderedParticles;
/**
* Storage of a potenital n-body decay
*/
class PrototypeVertex : public Base {
public:
/**
* Default Constructor
*/
PrototypeVertex() : npart(0), possibleOnShell(false) {}
/**
* Constructor
*/
PrototypeVertex(tPDPtr in, OrderedVertices out,
VertexBasePtr v, int n) :
incoming(in), outgoing(out), vertex(v), npart(n),
possibleOnShell(false) {}
/**
* Incoming particle
*/
tPDPtr incoming;
/**
* Outgoing particles
*/
OrderedVertices outgoing;
/**
* The vertex for the interaction
*/
VertexBasePtr vertex;
/**
* The parent of the vertex
*/
tPrototypeVertexPtr parent;
/**
* Number of particles
*/
unsigned int npart;
/**
* Outgoing particles
*/
mutable OrderedParticles outPart;
/**
* Can have on-shell intermediates
*/
bool possibleOnShell;
/**
* Increment the number of particles
*/
void incrementN(int in) {
npart += in;
if(parent) parent->incrementN(in);
}
/**
* Mass of the incoming particle
*/
Energy incomingMass() {
return incoming->mass();
}
/**
* Total mass of all the outgoing particles
*/
Energy outgoingMass() {
Energy mass(ZERO);
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
mass += it->second ?
it->second->outgoingMass() : it->first->mass();
}
return mass;
}
/**
* Total constituent mass of all the outgoing particles
*/
Energy outgoingConstituentMass() {
Energy mass(ZERO);
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
mass += it->second ?
it->second->outgoingConstituentMass() : it->first->constituentMass();
}
return mass;
}
/**
* Check the external particles
*/
bool checkExternal(bool first=true) {
if(outPart.empty()) setOutgoing();
if(first&&outPart.find(incoming)!=outPart.end()) return false;
bool output = true;
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
if(it->second&& !it->second->checkExternal(false)) output = false;
}
return output;
}
/**
* Set the outgoing particles
*/
void setOutgoing() const {
assert(outPart.empty());
for(OrderedVertices::const_iterator it = outgoing.begin();
it!=outgoing.end();++it) {
if(it->second) {
it->second->setOutgoing();
outPart.insert(it->second->outPart.begin(),
it->second->outPart.end());
}
else
outPart.insert(it->first);
}
}
/**
* Are there potential on-shell intermediates?
*/
bool canBeOnShell(unsigned int opt,Energy maxMass,bool first);
/**
* Check if same external particles
*/
bool sameDecay(const PrototypeVertex & x) const {
if(incoming != x.incoming) return false;
if(outPart.empty()) setOutgoing();
if(x.outPart.empty()) x.setOutgoing();
OrderedParticles::const_iterator cit = outPart.begin();
OrderedParticles::const_iterator cjt = x.outPart.begin();
if(x.npart!=npart) return false;
for(;cit!=outPart.end();++cit,++cjt) {
if(*cit!=*cjt) return false;
}
return true;
}
/**
* Create a \f$1\to2\f$ prototype
*/
static void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
std::queue<PrototypeVertexPtr> & prototypes);
/**
* Expand the prototypes by adding more legs
*/
static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
std::queue<PrototypeVertexPtr> & prototypes,
const set<PDPtr> & excluded);
/**
* Copy the whole structure with a new branching
*/
static PrototypeVertexPtr replicateTree(PrototypeVertexPtr parent,
PrototypeVertexPtr oldChild,
PrototypeVertexPtr & newChild);
};
/**
* Output to a stream
*/
inline ostream & operator<<(ostream & os, const PrototypeVertex & diag) {
os << diag.incoming->PDGName() << " -> ";
bool seq=false;
for(OrderedVertices::const_iterator it = diag.outgoing.begin();
it!=diag.outgoing.end();++it) {
os << it->first->PDGName() << " ";
if(it->second) seq = true;
}
os << " decays via "
<< diag.vertex->fullName() << " in a "
<< diag.npart << "-body decay\n";
if(!seq) return os;
os << "Followed by\n";
for(OrderedVertices::const_iterator it = diag.outgoing.begin();
it!=diag.outgoing.end();++it) {
if(it->second) os << *it->second;
}
return os;
}
/**
* Test whether two diagrams are identical.
*/
inline bool operator==(const PrototypeVertex & x, const PrototypeVertex & y) {
if(x.incoming != y.incoming) return false;
if(x.vertex != y.vertex) return false;
if(x.npart != y.npart) return false;
if(x.outgoing.size() != y.outgoing.size()) return false;
OrderedVertices::const_iterator xt = x.outgoing.begin();
OrderedVertices::const_iterator yt = y.outgoing.begin();
for(;xt!=x.outgoing.end();++xt,++yt) {
if(xt->first != yt->first) return false;
if(xt->second && yt->second) {
if(*(xt->second)==*(yt->second)) continue;
else return false;
}
else if(xt->second || yt->second)
return false;
}
return true;
}
/**
* A simple vertex for the N-body diagram
*/
struct NBVertex {
/**
* Constructor taking a prototype vertex as the arguments
*/
NBVertex(PrototypeVertexPtr proto = PrototypeVertexPtr() );
/**
* Incoming particle
*/
tPDPtr incoming;
/**
* Outgoing particles
*/
mutable OrderedParticles outgoing;
/**
* The vertices
*/
list<pair<PDPtr,NBVertex> > vertices;
/**
* The vertex
*/
VertexBasePtr vertex;
};
/**
* The NBDiagram struct contains information about a \f$1\ton\f$ decay
* that has been automatically generated.
*/
struct NBDiagram {
/**
- * Enumeration for the channel type
- */
- enum Channel {UNDEFINED = -1,
- TBchannel23=0, TBchannel13=1,
- TBchannel12=2, TBfourPoint=3};
-
- /**
- * Standard Constructor
- */
- NBDiagram() : channelType(UNDEFINED) {}
-
- /**
* Constructor taking a prototype vertex as the arguments*/
- NBDiagram(PrototypeVertexPtr proto);
+ NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());
/**
* Incoming particle
*/
tPDPtr incoming;
/**
* The type of channel
*/
- Channel channelType;
+ unsigned int channelType;
/**
* Outgoing particles
*/
mutable OrderedParticles outgoing;
/**
* The vertices
*/
list<pair<PDPtr,NBVertex> > vertices;
/**
* The vertex for the parent branching
*/
VertexBasePtr vertex;
/** Store colour flow at \f$N_c=3\f$ information */
mutable vector<CFPair> colourFlow;
/** Store colour flow at \f$N_c=\infty\f$ information */
mutable vector<CFPair> largeNcColourFlow;
};
-// /**
-// * Test whether this and x are the same decay
-// * @param x The other process to check
-// */
-// bool sameDecay(const TBDiagram & x) const {
-// if(ids[0] != x.ids[0]) return false;
-// bool used[4]={false,false,false,false};
-// for(unsigned int ix=1;ix<4;++ix) {
-// bool found=false;
-// for(unsigned int iy=1;iy<4;++iy) {
-// if(used[iy]) continue;
-// if(ids[ix]==x.ids[iy]) {
-// used[iy]=true;
-// found=true;
-// break;
-// }
-// }
-// if(!found) return false;
-// }
-// return true;
-// }
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
- * @param x The TBVertex
+ * @param x The NBVertex
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const NBVertex & x) {
os << x.incoming << x.outgoing << x.vertices << x.vertex;
return os;
}
/**
* Input operator to allow persistently written data to be read in
* @param is The input stream
- * @param x The TBVertex
+ * @param x The NBVertex
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
NBVertex & x) {
is >> x.incoming >> x.outgoing >> x.vertices >> x.vertex;
return is;
}
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
- * @param x The TBDiagram
+ * @param x The NBDiagram
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const NBDiagram & x) {
os << x.incoming << oenum(x.channelType) << x.outgoing << x.vertices << x.vertex
<< x.colourFlow << x.largeNcColourFlow;
return os;
}
/**
* Input operator to allow persistently written data to be read in
* @param is The input stream
- * @param x The TBDiagram
+ * @param x The NBDiagram
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
NBDiagram & x) {
is >> x.incoming >> ienum(x.channelType) >> x.outgoing >> x.vertices >> x.vertex
>> x.colourFlow >> x.largeNcColourFlow;
return is;
}
+/**
+ * Output a NBVertex to a stream
+ */
+inline ostream & operator<<(ostream & os, const NBVertex & vertex) {
+ os << vertex.incoming->PDGName() << " -> ";
+ bool seq=false;
+ for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
+ it!=vertex.vertices.end();++it) {
+ os << it->first->PDGName() << " ";
+ if(it->second.incoming) seq = true;
+ }
+ os << "via vertex " << vertex.vertex->fullName() << "\n";
+ if(!seq) return os;
+ os << "Followed by\n";
+ for(list<pair<PDPtr,NBVertex> >::const_iterator it=vertex.vertices.begin();
+ it!=vertex.vertices.end();++it) {
+ if(it->second.incoming) os << it->second;
+ }
+ return os;
+}
+
+/**
+ * Output a NBDiagram to a stream
+ */
+inline ostream & operator<<(ostream & os, const NBDiagram & diag) {
+ os << diag.incoming->PDGName() << " -> ";
+ for(OrderedParticles::const_iterator it=diag.outgoing.begin();
+ it!=diag.outgoing.end();++it) {
+ os << (**it).PDGName() << " ";
+ }
+ os << " has order " << diag.channelType << "\n";
+ os << "First decay " << diag.incoming->PDGName() << " -> ";
+ bool seq=false;
+ for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
+ it!=diag.vertices.end();++it) {
+ os << it->first->PDGName() << " ";
+ if(it->second.incoming) seq = true;
+ }
+ os << "via vertex " << diag.vertex->fullName() << "\n";
+ if(!seq) return os;
+ os << "Followed by\n";
+ for(list<pair<PDPtr,NBVertex> >::const_iterator it=diag.vertices.begin();
+ it!=diag.vertices.end();++it) {
+ if(it->second.incoming) os << it->second;
+ }
+ return os;
+}
+
}
#endif /* HERWIG_PrototypeVertex_H */

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:50 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805926
Default Alt Text
(99 KB)

Event Timeline