Page MenuHomeHEPForge
Files F7879826

No OneTemporary
Actions

View Options

This file is larger than 256 KB, so syntax highlighting was skipped.
diff --git a/Decay/General/GeneralTwoBodyDecayer.cc b/Decay/General/GeneralTwoBodyDecayer.cc
--- a/Decay/General/GeneralTwoBodyDecayer.cc
+++ b/Decay/General/GeneralTwoBodyDecayer.cc
@@ -1,1393 +1,1394 @@
// -*- C++ -*-
//
// GeneralTwoBodyDecayer.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the GeneralTwoBodyDecayer class.
//
#include "GeneralTwoBodyDecayer.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Utilities/Exception.h"
#include "Herwig++/Shower/Base/HardTree.h"
#include "Herwig++/Shower/Base/ShowerTree.h"
#include "Herwig++/Shower/Base/ShowerProgenitor.h"
#include "Herwig++/Shower/Base/ShowerParticle.h"
#include "Herwig++/Shower/Base/Branching.h"
using namespace Herwig;
ParticleVector GeneralTwoBodyDecayer::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;
}
void GeneralTwoBodyDecayer::doinit() {
DecayIntegrator::doinit();
assert( vertex_ );
assert( _incoming && _outgoing.size()==2);
vertex_->init();
//create phase space mode
tPDVector extpart(3);
extpart[0] = _incoming;
extpart[1] = _outgoing[0];
extpart[2] = _outgoing[1];
addMode(new_ptr(DecayPhaseSpaceMode(extpart, this)), _maxweight, vector<double>());
}
int GeneralTwoBodyDecayer::modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const {
long parentID = parent->id();
long id1 = children[0]->id();
long id2 = children[1]->id();
cc = false;
long out1 = _outgoing[0]->id();
long out2 = _outgoing[1]->id();
if( parentID == _incoming->id() &&
((id1 == out1 && id2 == out2) ||
(id1 == out2 && id2 == out1)) ) {
return 0;
}
else if(_incoming->CC() && parentID == _incoming->CC()->id()) {
cc = true;
if( _outgoing[0]->CC()) out1 = _outgoing[0]->CC()->id();
if( _outgoing[1]->CC()) out2 = _outgoing[1]->CC()->id();
if((id1 == out1 && id2 == out2) ||
(id1 == out2 && id2 == out1)) return 0;
}
return -1;
}
void GeneralTwoBodyDecayer::
colourConnections(const Particle & parent,
const ParticleVector & out) const {
PDT::Colour incColour(parent.data().iColour());
PDT::Colour outaColour(out[0]->data().iColour());
PDT::Colour outbColour(out[1]->data().iColour());
//incoming colour singlet
if(incColour == PDT::Colour0) {
// colour triplet-colourantitriplet
if((outaColour == PDT::Colour3 && outbColour == PDT::Colour3bar) ||
(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3)) {
bool ac(out[0]->id() < 0);
out[0]->colourNeighbour(out[1],!ac);
}
//colour octet
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour8) {
out[0]->colourNeighbour(out[1]);
out[0]->antiColourNeighbour(out[1]);
}
// colour singlets
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour0) {
}
// unknown
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour singlet in "
<< "GeneralTwoBodyDecayer::colourConnections "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
//incoming colour triplet
else if(incColour == PDT::Colour3) {
// colour triplet + singlet
if(outaColour == PDT::Colour3 && outbColour == PDT::Colour0) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
//opposite order
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
// octet + triplet
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->antiColourNeighbour(out[0]);
}
//opposite order
else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour8) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[0]->antiColourNeighbour(out[1]);
}
else if(outaColour == PDT::Colour3bar && outaColour == PDT::Colour3bar) {
tColinePtr col[2] = {ColourLine::create(out[0],true),
ColourLine::create(out[1],true)};
parent.colourLine()->setSinkNeighbours(col[0],col[1]);
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour triplet in "
<< "GeneralTwoBodyDecayer::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
// incoming colour anti triplet
else if(incColour == PDT::Colour3bar) {
// colour antitriplet +singlet
if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour0) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
//opposite order
else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3bar) {
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
//octet + antitriplet
else if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour8) {
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[0]->colourNeighbour(out[1]);
}
//opposite order
else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3bar) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[1]->colourNeighbour(out[0]);
}
else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
tColinePtr col[2] = {ColourLine::create(out[0]),
ColourLine::create(out[1])};
parent.antiColourLine()->setSourceNeighbours(col[0],col[1]);
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour antitriplet "
<< "in GeneralTwoBodyDecayer::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
//incoming colour octet
else if(incColour == PDT::Colour8) {
// triplet-antitriplet
if(outaColour == PDT::Colour3&&outbColour == PDT::Colour3bar) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
// opposite order
else if(outbColour == PDT::Colour3&&outaColour == PDT::Colour3bar) {
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
}
// neutral octet
else if(outaColour == PDT::Colour0&&outbColour == PDT::Colour8) {
out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
else if(outbColour == PDT::Colour0&&outaColour == PDT::Colour8) {
out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour octet "
<< "in GeneralTwoBodyDecayer::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else if(incColour == PDT::Colour6) {
if(outaColour == PDT::Colour3 && outbColour == PDT::Colour3) {
tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(tempParent->colourInfo());
tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[0]);
line1->addColoured(dynamic_ptr_cast<tPPtr>(out[0]));
tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->colourLines()[1]);
line2->addColoured(dynamic_ptr_cast<tPPtr>(out[1]));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour sextet "
<< "in GeneralTwoBodyDecayer::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else if(incColour == PDT::Colour6bar) {
if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3bar) {
tPPtr tempParent = const_ptr_cast<tPPtr>(&parent);
Ptr<MultiColour>::pointer parentColour =
dynamic_ptr_cast<Ptr<MultiColour>::pointer>
(tempParent->colourInfo());
tColinePtr line1 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[0]);
line1->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[0]));
tColinePtr line2 = const_ptr_cast<tColinePtr>(parentColour->antiColourLines()[1]);
line2->addAntiColoured(dynamic_ptr_cast<tPPtr>(out[1]));
}
else
throw Exception() << "Unknown outgoing colours for decaying "
<< "colour anti-sextet "
<< "in GeneralTwoBodyDecayer::colourConnections() "
<< outaColour << " " << outbColour
<< Exception::runerror;
}
else
throw Exception() << "Unknown incoming colour in "
<< "GeneralTwoBodyDecayer::colourConnections() "
<< incColour
<< Exception::runerror;
}
bool GeneralTwoBodyDecayer::twoBodyMEcode(const DecayMode & dm, int & mecode,
double & coupling) const {
assert(dm.parent()->id() == _incoming->id());
ParticleMSet::const_iterator pit = dm.products().begin();
long id1 = (*pit)->id();
++pit;
long id2 = (*pit)->id();
long id1t(_outgoing[0]->id()), id2t(_outgoing[1]->id());
mecode = -1;
coupling = 1.;
if( id1 == id1t && id2 == id2t ) {
return true;
}
else if( id1 == id2t && id2 == id1t ) {
return false;
}
else
assert(false);
return false;
}
void GeneralTwoBodyDecayer::persistentOutput(PersistentOStream & os) const {
os << vertex_ << _incoming << _outgoing << _maxweight << ounit(pTmin_,GeV)
<< coupling_ << incomingVertex_ << outgoingVertices_ << fourPointVertex_;
}
void GeneralTwoBodyDecayer::persistentInput(PersistentIStream & is, int) {
is >> vertex_ >> _incoming >> _outgoing >> _maxweight >> iunit(pTmin_,GeV)
>> coupling_ >> incomingVertex_ >> outgoingVertices_ >> fourPointVertex_;
}
AbstractClassDescription<GeneralTwoBodyDecayer>
GeneralTwoBodyDecayer::initGeneralTwoBodyDecayer;
// Definition of the static class description member.
void GeneralTwoBodyDecayer::Init() {
static ClassDocumentation<GeneralTwoBodyDecayer> documentation
("This class is designed to be a base class for all 2 body decays"
"in a general model");
static Parameter<GeneralTwoBodyDecayer,Energy> interfacepTmin
("pTmin",
"Minimum transverse momentum from gluon radiation",
&GeneralTwoBodyDecayer::pTmin_, GeV, 1.0*GeV, 0.0*GeV, 10.0*GeV,
false, false, Interface::limited);
static Reference<GeneralTwoBodyDecayer,ShowerAlpha> interfaceCoupling
("Coupling",
"Object for the coupling in the generation of hard radiation",
&GeneralTwoBodyDecayer::coupling_, false, false, true, false, false);
}
double GeneralTwoBodyDecayer::brat(const DecayMode &, const Particle & p,
double oldbrat) const {
ParticleVector children = p.children();
if( children.size() != 2 || !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()) );
Energy width = p.data().widthGenerator()->width(p.data(), scale);
return pwidth/width;
}
void GeneralTwoBodyDecayer::doinitrun() {
+ vertex_->initrun();
DecayIntegrator::doinitrun();
for(unsigned int ix=0;ix<numberModes();++ix) {
double fact = pow(1.5,int(mode(ix)->externalParticles(0)->iSpin())-1);
mode(ix)->setMaxWeight(fact*mode(ix)->maxWeight());
}
}
double GeneralTwoBodyDecayer::colourFactor(tcPDPtr in, tcPDPtr out1,
tcPDPtr out2) const {
// identical particle symmetry factor
double output = out1->id()==out2->id() ? 0.5 : 1.;
// colour neutral incoming particle
if(in->iColour()==PDT::Colour0) {
// both colour neutral
if(out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour0)
output *= 1.;
// colour triplet/ antitriplet
else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
output *= 3.;
}
// colour octet colour octet
else if(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour8 ) {
output *= 8.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour neutral particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
// triplet
else if(in->iColour()==PDT::Colour3) {
// colour triplet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3) ||
(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour0) ) {
output *= 1.;
}
// colour triplet + octet
else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3) ||
(out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour8) ) {
output *= 4./3.;
}
// colour anti triplet anti triplet
else if(out1->iColour()==PDT::Colour3bar &&
out2->iColour()==PDT::Colour3bar) {
output *= 2.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour triplet particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
// anti triplet
else if(in->iColour()==PDT::Colour3bar) {
// colour anti triplet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour3bar ) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour0 ) ) {
output *= 1.;
}
// colour anti triplet + octet
else if((out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour3bar ) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour8 ) ) {
output *= 4./3.;
}
// colour triplet triplet
else if(out1->iColour()==PDT::Colour3 &&
out2->iColour()==PDT::Colour3) {
output *= 2.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti triplet particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour8) {
// colour octet + neutral
if((out1->iColour()==PDT::Colour0 && out2->iColour()==PDT::Colour8 ) ||
(out1->iColour()==PDT::Colour8 && out2->iColour()==PDT::Colour0 ) ) {
output *= 1.;
}
// colour triplet/antitriplet
else if((out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3bar) ||
(out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3 ) ) {
output *= 0.5;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour octet particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour6) {
// colour sextet -> triplet triplet
if( out1->iColour()==PDT::Colour3 && out2->iColour()==PDT::Colour3 ) {
output *= 1.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour sextet particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else if(in->iColour()==PDT::Colour6bar) {
// colour sextet -> triplet triplet
if( out1->iColour()==PDT::Colour3bar && out2->iColour()==PDT::Colour3bar ) {
output *= 1.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti-sextet particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
}
else
throw Exception() << "Unknown colour "
<< in->iColour() << " for the decaying particle in "
<< "GeneralTwoBodyDecayer::colourFactor() for "
<< in->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< Exception::runerror;
return output;
}
Energy GeneralTwoBodyDecayer::partialWidth(PMPair inpart, PMPair outa,
PMPair outb) const {
// select the number of the mode
tPDVector children;
children.push_back(const_ptr_cast<PDPtr>(outa.first));
children.push_back(const_ptr_cast<PDPtr>(outb.first));
bool cc;
int nmode=modeNumber(cc,inpart.first,children);
tcPDPtr newchild[2] = {mode(nmode)->externalParticles(1),
mode(nmode)->externalParticles(2)};
// make the particles
Lorentz5Momentum pparent = Lorentz5Momentum(inpart.second);
PPtr parent = inpart.first->produceParticle(pparent);
Lorentz5Momentum pout[2];
double ctheta,phi;
Kinematics::generateAngles(ctheta,phi);
Kinematics::twoBodyDecay(pparent, outa.second, outb.second,
ctheta, phi,pout[0],pout[1]);
if( ( !cc && outa.first!=newchild[0]) ||
( cc && !(( outa.first->CC() && outa.first->CC() == newchild[0])||
( !outa.first->CC() && outa.first == newchild[0]) )))
swap(pout[0],pout[1]);
ParticleVector decay;
decay.push_back(newchild[0]->produceParticle(pout[0]));
decay.push_back(newchild[1]->produceParticle(pout[1]));
double me = me2(-1,*parent,decay,Initialize);
Energy pcm = Kinematics::pstarTwoBodyDecay(inpart.second,
outa.second, outb.second);
return me/(8.*Constants::pi)*pcm;
}
void GeneralTwoBodyDecayer::setDecayInfo(PDPtr incoming,PDPair outgoing,
VertexBasePtr vertex, VertexBasePtr inV,
const vector<VertexBasePtr> & outV,
VertexBasePtr fourV) {
_incoming=incoming;
_outgoing.clear();
_outgoing.push_back(outgoing.first );
_outgoing.push_back(outgoing.second);
vertex_ = vertex;
incomingVertex_ = inV;
outgoingVertices_ = outV;
fourPointVertex_ = fourV;
}
HardTreePtr GeneralTwoBodyDecayer::generateHardest(ShowerTreePtr tree) {
// ignore effective vertices
if (vertex_ && (vertex_->orderInGem()+vertex_->orderInGs())>1)
return HardTreePtr();
// search for coloured particles
bool colouredParticles=false;
vector<ShowerProgenitorPtr> Progenitors = tree->extractProgenitors();
for (unsigned int it=0; it<Progenitors.size(); ++it){
if (Progenitors[it]->progenitor()->dataPtr()->coloured()){
colouredParticles=true;
break;
}
}
// if no coloured particles return
if ( !colouredParticles ) return HardTreePtr();
// check exactly two outgoing particles
if (tree->outgoingLines().size()!=2)
assert(false);
// for decay b -> a c
// set progenitors
ShowerProgenitorPtr
cProgenitor = tree->outgoingLines(). begin()->first,
aProgenitor = tree->outgoingLines().rbegin()->first;
// get the decaying particle
ShowerProgenitorPtr bProgenitor = tree->incomingLines().begin()->first;
// identify which dipoles are required
vector<dipoleType> dipoles;
if(!identifyDipoles(dipoles,aProgenitor,bProgenitor,cProgenitor)) {
return HardTreePtr();
}
Energy trialpT = pTmin_;
LorentzRotation eventFrame;
vector<Lorentz5Momentum> momenta;
vector<Lorentz5Momentum> trialMomenta(4);
ShowerProgenitorPtr finalEmitter, finalSpectator;
ShowerProgenitorPtr trialEmitter, trialSpectator;
for (int i=0; i<int(dipoles.size()); ++i){
// assign emitter and spectator based on current dipole
if (dipoles[i]==FFc || dipoles[i]==IFc || dipoles[i]==IFbc){
trialEmitter = cProgenitor;
trialSpectator = aProgenitor;
}
else if (dipoles[i]==FFa || dipoles[i]==IFa || dipoles[i]==IFba){
trialEmitter = aProgenitor;
trialSpectator = cProgenitor;
}
// find rotation from lab to frame with the spectator along -z
LorentzRotation trialEventFrame = ( bProgenitor->progenitor()->momentum()
.findBoostToCM() );
Lorentz5Momentum pspectator = (trialEventFrame*
trialSpectator->progenitor()->momentum());
trialEventFrame.rotateZ( -pspectator.phi() );
trialEventFrame.rotateY( -pspectator.theta() - Constants::pi );
// invert it
trialEventFrame.invert();
// try to generate an emission
pT_ = pTmin_;
vector<Lorentz5Momentum> trialMomenta
= hardMomenta(bProgenitor, trialEmitter, trialSpectator, dipoles, i);
// select dipole which gives highest pT emission
if(pT_>trialpT){
trialpT = pT_;
momenta = trialMomenta;
eventFrame = trialEventFrame;
finalEmitter = trialEmitter;
finalSpectator = trialSpectator;
}
}
pT_ = trialpT;
// if no emission return
if(momenta.empty()) {
bProgenitor->maximumpT(pTmin_,ShowerInteraction::QCD);
aProgenitor->maximumpT(pTmin_,ShowerInteraction::QCD);
cProgenitor->maximumpT(pTmin_,ShowerInteraction::QCD);
return HardTreePtr();
}
// rotate momenta back to the lab
for(unsigned int ix=0;ix<momenta.size();++ix) {
momenta[ix] *= eventFrame;
}
// set maximum pT for subsequent branchings
bProgenitor ->maximumpT(pT_,ShowerInteraction::QCD);
finalEmitter ->maximumpT(pT_,ShowerInteraction::QCD);
finalSpectator->maximumpT(pT_,ShowerInteraction::QCD);
// get ParticleData objects
tcPDPtr b = bProgenitor ->progenitor()->dataPtr();
tcPDPtr e = finalEmitter ->progenitor()->dataPtr();
tcPDPtr s = finalSpectator->progenitor()->dataPtr();
tcPDPtr gluon = getParticleData(ParticleID::g);
// create new ShowerParticles
ShowerParticlePtr emitter (new_ptr(ShowerParticle(e, true )));
ShowerParticlePtr spectator(new_ptr(ShowerParticle(s, true )));
ShowerParticlePtr gauge (new_ptr(ShowerParticle(gluon, true )));
ShowerParticlePtr incoming (new_ptr(ShowerParticle(b, false)));
ShowerParticlePtr parent (new_ptr(ShowerParticle(e, true )));
// set momenta
emitter ->set5Momentum(momenta[1]);
spectator->set5Momentum(momenta[2]);
gauge ->set5Momentum(momenta[3]);
incoming ->set5Momentum(bProgenitor->progenitor()->momentum());
Lorentz5Momentum parentMomentum(momenta[1]+momenta[3]);
parentMomentum.rescaleMass();
parent->set5Momentum(parentMomentum);
// create the vectors of HardBranchings to create the HardTree:
vector<HardBranchingPtr> spaceBranchings,allBranchings;
// incoming particle b
spaceBranchings.push_back(new_ptr(HardBranching(incoming,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Incoming)));
// outgoing particles from hard emission:
HardBranchingPtr spectatorBranch(new_ptr(HardBranching(spectator,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
HardBranchingPtr emitterBranch(new_ptr(HardBranching(parent,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
emitterBranch->addChild(new_ptr(HardBranching(emitter,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
emitterBranch->addChild(new_ptr(HardBranching(gauge,SudakovPtr(),
HardBranchingPtr(),
HardBranching::Outgoing)));
if (emitterBranch->branchingParticle()->dataPtr()->hasColour()
&& (! emitterBranch->branchingParticle()->dataPtr()->hasAntiColour())) {
emitterBranch->type(ShowerPartnerType::QCDColourLine);
}
else if (emitterBranch->branchingParticle()->dataPtr()->hasAntiColour()
&& (! emitterBranch->branchingParticle()->dataPtr()->hasColour())) {
emitterBranch->type(ShowerPartnerType::QCDAntiColourLine);
}
else
emitterBranch->type(UseRandom::rndbool() ?
ShowerPartnerType::QCDColourLine :
ShowerPartnerType::QCDAntiColourLine);
allBranchings.push_back(spaceBranchings[0]);
allBranchings.push_back(emitterBranch);
allBranchings.push_back(spectatorBranch);
// make the HardTree from the HardBranching vectors.
HardTreePtr hardtree = new_ptr(HardTree(allBranchings,spaceBranchings,
ShowerInteraction::QCD));
// connect the particles with the branchings in the HardTree
hardtree->connect( bProgenitor ->progenitor(), spaceBranchings[0] );
hardtree->connect( finalEmitter ->progenitor(), allBranchings[1] );
hardtree->connect( finalSpectator->progenitor(), allBranchings[2] );
// set up colour lines
vector<ColinePtr> newline;
getColourLines(newline, hardtree, bProgenitor);
// return the tree
return hardtree;
}
double GeneralTwoBodyDecayer::threeBodyME(const int , const Particle &,
const ParticleVector &,MEOption) {
throw Exception() << "Base class GeneralTwoBodyDecayer::threeBodyME() "
<< "called, should have an implementation in the inheriting class"
<< Exception::runerror;
return 0.;
}
vector<Lorentz5Momentum> GeneralTwoBodyDecayer::hardMomenta(const ShowerProgenitorPtr &in,
const ShowerProgenitorPtr &emitter,
const ShowerProgenitorPtr &spectator,
const vector<dipoleType> &dipoles,
int i) {
double C = 6.3;
double ymax = 10.;
double ymin = -ymax;
// get masses of the particles
mb_ = in ->progenitor()->momentum().mass();
e_ = emitter ->progenitor()->momentum().mass()/mb_;
s_ = spectator->progenitor()->momentum().mass()/mb_;
e2_ = sqr(e_);
s2_ = sqr(s_);
vector<Lorentz5Momentum> particleMomenta (4);
Energy2 lambda = sqr(mb_)*sqrt(1.+sqr(s2_)+sqr(e2_)-2.*s2_-2.*e2_-2.*s2_*e2_);
// calculate A
double A = (ymax-ymin)*C*(coupling()->overestimateValue()/(2.*Constants::pi));
Energy pTmax = mb_*(sqr(1.-s_)-e2_)/(2.*(1.-s_));
// if no possible branching return
if ( pTmax < pTmin_ ) {
particleMomenta.clear();
return particleMomenta;
}
while (pTmax >= pTmin_) {
// generate pT, y and phi values
Energy pT = pTmax*pow(UseRandom::rnd(),(1./A));
if (pT < pTmin_) {
particleMomenta.clear();
break;
}
double phi = UseRandom::rnd()*Constants::twopi;
double y = ymin+UseRandom::rnd()*(ymax-ymin);
double weight[2] = {0.,0.};
double xs[2], xe[2], xe_z[2], xg;
for (unsigned int j=0; j<2; j++) {
// check if the momenta are physical
if (!calcMomenta(j, pT, y, phi, xg, xs[j], xe[j], xe_z[j], particleMomenta))
continue;
// check if point lies within phase space
if (!psCheck(xg, xs[j]))
continue;
// decay products for 3 body decay
PPtr inpart = in ->progenitor()->dataPtr()->produceParticle(particleMomenta[0]);
ParticleVector decay3;
decay3.push_back(emitter ->progenitor()->dataPtr()->produceParticle(particleMomenta[1]));
decay3.push_back(spectator->progenitor()->dataPtr()->produceParticle(particleMomenta[2]));
decay3.push_back(getParticleData(ParticleID::g )->produceParticle(particleMomenta[3]));
// decay products for 2 body decay
Lorentz5Momentum p1(ZERO,ZERO, lambda/2./mb_,(mb_/2.)*(1.+e2_-s2_),mb_*e_);
Lorentz5Momentum p2(ZERO,ZERO,-lambda/2./mb_,(mb_/2.)*(1.+s2_-e2_),mb_*s_);
ParticleVector decay2;
decay2.push_back(emitter ->progenitor()->dataPtr()->produceParticle(p1));
decay2.push_back(spectator->progenitor()->dataPtr()->produceParticle(p2));
if (decay2[0]->dataPtr()->iSpin()!=PDT::Spin1Half &&
decay2[1]->dataPtr()->iSpin()==PDT::Spin1Half) swap(decay2[0], decay2[1]);
// calculate matrix element ratio R/B
double meRatio = matrixElementRatio(*inpart,decay2,decay3,Initialize);
// calculate dipole factor
InvEnergy2 dipoleSum = ZERO;
InvEnergy2 numerator = ZERO;
for (int k=0; k<int(dipoles.size()); ++k){
InvEnergy2 dipole = abs(calculateDipole(dipoles[k],*inpart,decay3,dipoles[i]));
dipoleSum += dipole;
if (k==i) numerator = dipole;
}
meRatio *= numerator/dipoleSum;
// calculate jacobian
Energy2 denom = (mb_-particleMomenta[3].e())*particleMomenta[2].vect().mag() -
particleMomenta[2].e()*particleMomenta[3].z();
InvEnergy2 J = (particleMomenta[2].vect().mag2())/(lambda*denom);
// calculate weight
weight[j] = meRatio*fabs(sqr(pT)*J)*coupling()->ratio(pT*pT)/C/Constants::twopi;
}
// accept point if weight > R
if (weight[0] + weight[1] > UseRandom::rnd()) {
if (weight[0] > (weight[0] + weight[1])*UseRandom::rnd()) {
particleMomenta[1].setE( (mb_/2.)*xe [0]);
particleMomenta[1].setZ( (mb_/2.)*xe_z[0]);
particleMomenta[2].setE( (mb_/2.)*xs [0]);
particleMomenta[2].setZ(-(mb_/2.)*sqrt(sqr(xs[0])-4.*s2_));
}
else {
particleMomenta[1].setE( (mb_/2.)*xe [1]);
particleMomenta[1].setZ( (mb_/2.)*xe_z[1]);
particleMomenta[2].setE( (mb_/2.)*xs [1]);
particleMomenta[2].setZ(-(mb_/2.)*sqrt(sqr(xs[1])-4.*s2_));
}
pT_ = pT;
break;
}
// if there's no branching lower the pT
pTmax = pT;
}
return particleMomenta;
}
double GeneralTwoBodyDecayer::matrixElementRatio(const Particle & inpart,
const ParticleVector & decay2,
const ParticleVector & decay3,
MEOption meopt) {
// calculate R/B
double B = me2 (0, inpart, decay2, meopt);
double R = threeBodyME(0, inpart, decay3, meopt);
return R/B;
}
bool GeneralTwoBodyDecayer::calcMomenta(int j, Energy pT, double y, double phi,
double& xg, double& xs, double& xe, double& xe_z,
vector<Lorentz5Momentum>& particleMomenta){
// calculate xg
xg = 2.*pT*cosh(y) / mb_;
if (xg>(1. - sqr(e_ + s_)) || xg<0.) return false;
// calculate the two values of xs
double xT = 2.*pT / mb_;
double A = 4.-4.*xg+sqr(xT);
double B = 4.*(3.*xg-2.+2.*e2_-2.*s2_-sqr(xg)-xg*e2_+xg*s2_);
double L = 1.+sqr(s2_)+sqr(e2_)-2.*s2_-2.*e2_-2.*s2_*e2_;
double det = 16.*( -L*sqr(xT)+pow(xT,4)*s2_+2.*xg*sqr(xT)*(1.-s2_-e2_)+
L*sqr(xg)-sqr(xg*xT)*(1.+s2_)+pow(xg,4)+
2.*pow(xg,3)*(-1.+s2_+e2_) );
if (det<0.) return false;
if (j==0) xs = (-B+sqrt(det))/(2.*A);
if (j==1) xs = (-B-sqrt(det))/(2.*A);
// check value of xs is physical
if (xs>(1.+s2_-e2_) || xs<2.*s_) return false;
// calculate xe
xe = 2.-xs-xg;
// check value of xe is physical
if (xe>(1.+e2_-s2_) || xe<2.*e_) return false;
// calculate xe_z
double epsilon_p = -sqrt(sqr(xs)-4.*s2_)+xT*sinh(y)+sqrt(sqr(xe)-4.*e2_-sqr(xT));
double epsilon_m = -sqrt(sqr(xs)-4.*s2_)+xT*sinh(y)-sqrt(sqr(xe)-4.*e2_-sqr(xT));
// find direction of emitter
if (fabs(epsilon_p) < 1.e-10) xe_z = sqrt(sqr(xe)-4.*e2_-sqr(xT));
else if (fabs(epsilon_m) < 1.e-10) xe_z = -sqrt(sqr(xe)-4.*e2_-sqr(xT));
else return false;
// check the emitter is on shell
if (fabs((sqr(xe)-sqr(xT)-sqr(xe_z)-4.*e2_))>1.e-10) return false;
// calculate 4 momenta
particleMomenta[0].setE ( mb_);
particleMomenta[0].setX ( ZERO);
particleMomenta[0].setY ( ZERO);
particleMomenta[0].setZ ( ZERO);
particleMomenta[0].setMass( mb_);
particleMomenta[1].setE ( mb_*xe/2.);
particleMomenta[1].setX (-pT*cos(phi));
particleMomenta[1].setY (-pT*sin(phi));
particleMomenta[1].setZ ( mb_*xe_z/2.);
particleMomenta[1].setMass( mb_*e_);
particleMomenta[2].setE ( mb_*xs/2.);
particleMomenta[2].setX ( ZERO);
particleMomenta[2].setY ( ZERO);
particleMomenta[2].setZ (-mb_*sqrt(sqr(xs)-4.*s2_)/2.);
particleMomenta[2].setMass( mb_*s_);
particleMomenta[3].setE ( pT*cosh(y));
particleMomenta[3].setX ( pT*cos(phi));
particleMomenta[3].setY ( pT*sin(phi));
particleMomenta[3].setZ ( pT*sinh(y));
particleMomenta[3].setMass( ZERO);
return true;
}
bool GeneralTwoBodyDecayer::psCheck(const double xg, const double xs) {
// check is point is in allowed region of phase space
double xe_star = (1.-s2_+e2_-xg)/sqrt(1.-xg);
double xg_star = xg/sqrt(1.-xg);
if ((sqr(xe_star)-4.*e2_) < 1e-10) return false;
double xs_max = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)+xg_star))/ 4.;
double xs_min = (4.+4.*s2_-sqr(xe_star+xg_star)+
sqr(sqrt(sqr(xe_star)-4.*e2_)-xg_star))/ 4.;
if (xs < xs_min || xs > xs_max) return false;
return true;
}
double GeneralTwoBodyDecayer::colourCoeff(const PDT::Colour emitter,
const PDT::Colour spectator,
const PDT::Colour other){
// calculate the colour factor of the dipole
double numerator=1.;
double denominator=1.;
if (emitter!=PDT::Colour0 && spectator!=PDT::Colour0 && other!=PDT::Colour0){
if (emitter ==PDT::Colour3 || emitter ==PDT::Colour3bar) numerator=-4./3;
else if (emitter ==PDT::Colour8) numerator=-3. ;
denominator=-1.*numerator;
if (spectator==PDT::Colour3 || spectator==PDT::Colour3bar) numerator-=4./3;
else if (spectator==PDT::Colour8) numerator-=3. ;
if (other ==PDT::Colour3 || other ==PDT::Colour3bar) numerator+=4./3;
else if (other ==PDT::Colour8) numerator+=3. ;
numerator*=(-1./2.);
}
if (emitter==PDT::Colour3 || emitter== PDT::Colour3bar) numerator*=4./3.;
else if (emitter==PDT::Colour8 && spectator!=PDT::Colour8) numerator*=3.;
else if (emitter==PDT::Colour8 && spectator==PDT::Colour8) numerator*=6.;
return (numerator/denominator);
}
InvEnergy2 GeneralTwoBodyDecayer::calculateDipole(const dipoleType & dipoleId,
const Particle & inpart,
const ParticleVector & decay3,
const dipoleType & emittingDipole){
// calculate dipole for decay b->ac
InvEnergy2 dipole = ZERO;
double xe = 2.*decay3[0]->momentum().e()/mb_;
double xs = 2.*decay3[1]->momentum().e()/mb_;
double xg = 2.*decay3[2]->momentum().e()/mb_;
double coeff = 8.*Constants::pi*coupling()->value(mb_*mb_);
// radiation from b with initial-final connection
if (dipoleId==IFba || dipoleId==IFbc){
dipole = -2./sqr(mb_*xg);
dipole *= colourCoeff(inpart.dataPtr()->iColour(), decay3[0]->dataPtr()->iColour(),
decay3[1]->dataPtr()->iColour());
}
// radiation from a/c with initial-final connection
else if ((dipoleId==IFa &&
(emittingDipole==IFba || emittingDipole==IFa || emittingDipole==FFa)) ||
(dipoleId==IFc &&
(emittingDipole==IFbc || emittingDipole==IFc || emittingDipole==FFc))){
double z = 1. - xg/(1.-s2_+e2_);
dipole = (-2.*e2_/sqr(1.-xs+s2_-e2_)/sqr(mb_) + (1./(1.-xs+s2_-e2_)/sqr(mb_))*
(2./(1.-z)-dipoleSpinFactor(decay3[0],z)));
dipole *= colourCoeff(decay3[0]->dataPtr()->iColour(),inpart.dataPtr()->iColour(),
decay3[1]->dataPtr()->iColour());
}
else if (dipoleId==IFa || dipoleId==IFc){
double z = 1. - xg/(1.-e2_+s2_);
dipole = (-2.*s2_/sqr(1.-xe+e2_-s2_)/sqr(mb_)+(1./(1.-xe+e2_-s2_)/sqr(mb_))*
(2./(1.-z)-dipoleSpinFactor(decay3[1],z)));
dipole *= colourCoeff(decay3[1]->dataPtr()->iColour(),inpart.dataPtr()->iColour(),
decay3[0]->dataPtr()->iColour());
}
// radiation from a/c with final-final connection
else if ((dipoleId==FFa &&
(emittingDipole==IFba || emittingDipole==IFa || emittingDipole==FFa)) ||
(dipoleId==FFc &&
(emittingDipole==IFbc || emittingDipole==IFc || emittingDipole==FFc))){
double z = 1. + ((xe-1.+s2_-e2_)/(xs-2.*s2_));
dipole = (1./(1.-xs+s2_-e2_)/sqr(mb_))*((2./(1.-z))-dipoleSpinFactor(decay3[0],z)-
(2.*e2_/(1.+s2_-e2_-xs)) );
dipole *= colourCoeff(decay3[0]->dataPtr()->iColour(),
decay3[1]->dataPtr()->iColour(),
inpart.dataPtr()->iColour());
}
else if (dipoleId==FFa || dipoleId==FFc) {
double z = 1. + ((xs-1.+e2_-s2_)/(xe-2.*e2_));
dipole = (1./(1.-xe+e2_-s2_)/sqr(mb_))*((2./(1.-z))-dipoleSpinFactor(decay3[1],z)-
(2.*s2_/(1.+e2_-s2_-xe)) );
dipole *= colourCoeff(decay3[1]->dataPtr()->iColour(),
decay3[0]->dataPtr()->iColour(),
inpart.dataPtr()->iColour());
}
dipole *= coeff;
return dipole;
}
double GeneralTwoBodyDecayer::dipoleSpinFactor(const PPtr & emitter, double z){
// calculate the spin dependent component of the dipole
if (emitter->dataPtr()->iSpin()==PDT::Spin0)
return 2.;
else if (emitter->dataPtr()->iSpin()==PDT::Spin1Half)
return (1. + z);
else if (emitter->dataPtr()->iSpin()==PDT::Spin1)
return (2.*z*(1.-z) - 1./(1.-z) + 1./z -2.);
return 0.;
}
const vector<DVector> & GeneralTwoBodyDecayer::getColourFactors(const Particle & inpart,
const ParticleVector & decay,
unsigned int & nflow){
// calculate the colour factors for the three-body decay
vector<int> sing,trip,atrip,oct;
for(unsigned int it=0;it<decay.size();++it) {
if (decay[it]->dataPtr()->iColour() == PDT::Colour0 ) sing. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour3 ) trip. push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour3bar ) atrip.push_back(it);
else if(decay[it]->dataPtr()->iColour() == PDT::Colour8 ) oct. push_back(it);
}
// require at least one gluon
assert(oct.size()>=1);
// identical particle symmetry factor
double symFactor=1.;
if (( sing.size()==2 && decay[ sing[0]]->id()==decay[ sing[1]]->id()) ||
( trip.size()==2 && decay[ trip[0]]->id()==decay[ trip[1]]->id()) ||
(atrip.size()==2 && decay[atrip[0]]->id()==decay[atrip[1]]->id()) ||
( oct.size()==2 && decay[ oct[0]]->id()==decay[ oct[1]]->id()))
symFactor/=2.;
else if (oct.size()==3 &&
decay[oct[0]]->id()==decay[oct[1]]->id() &&
decay[oct[0]]->id()==decay[oct[2]]->id())
symFactor/=6.;
colour_ = vector<DVector>(1,DVector(1,symFactor*1.));
// decaying colour singlet
if(inpart.dataPtr()->iColour() == PDT::Colour0) {
if(trip.size()==1 && atrip.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4.));
}
else if (oct.size()==3){
nflow = 1.;
colour_ = vector<DVector>(1,DVector(1,symFactor*24.));
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour scalar particle in "
<< "GeneralTwoBodyDecayer::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3) {
if(trip.size()==1 && sing.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
}
else if(trip.size()==1 && oct.size()==2) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = symFactor*16./9.; colour_[0][1] = -symFactor*2./9.;
colour_[1][0] = -symFactor*2./9.; colour_[1][1] = symFactor*16./9.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour triplet particle in "
<< "GeneralTwoBodyDecayer::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour anti-triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3bar) {
if(atrip.size()==1 && sing.size()==1 && oct.size()==1) {
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*4./3.));
}
else if(atrip.size()==1 && oct.size()==2){
nflow = 2;
colour_.clear();
colour_ .resize(2,DVector(2,0.));
colour_[0][0] = symFactor*16./9.; colour_[0][1] = -symFactor*2./9.;
colour_[1][0] = -symFactor*2./9.; colour_[1][1] = symFactor*16./9.;
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour anti-triplet particle in "
<< "GeneralTwoBodyDecayer::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
// decaying colour octet
else if(inpart.dataPtr()->iColour() == PDT::Colour8) {
if(oct.size()==1 && trip.size()==1 && atrip.size()==1) {
nflow = 2;
colour_.clear();
colour_.resize(2,DVector(2,0.));
colour_[0][0] = symFactor*2./3. ; colour_[0][1] = -symFactor*1./12.;
colour_[1][0] = -symFactor*1./12.; colour_[1][1] = symFactor*2./3. ;
}
else if (oct.size()==2 && sing.size()==1){
nflow = 1;
colour_ = vector<DVector>(1,DVector(1,symFactor*3.));
}
else
throw Exception() << "Unknown colour for the outgoing particles"
<< " for decay colour octet particle in "
<< "GeneralTwoBodyDecayer::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
}
else
throw Exception() << "Unknown colour for the decaying particle in "
<< "GeneralTwoBodyDecayer::getColourFactors() for "
<< inpart. dataPtr()->PDGName() << " -> "
<< decay[0]->dataPtr()->PDGName() << " "
<< decay[1]->dataPtr()->PDGName() << " "
<< decay[2]->dataPtr()->PDGName()
<< Exception::runerror;
return colour_;
}
bool GeneralTwoBodyDecayer::identifyDipoles(vector<dipoleType> & dipoles,
ShowerProgenitorPtr & aProgenitor,
ShowerProgenitorPtr & bProgenitor,
ShowerProgenitorPtr & cProgenitor) const {
PDT::Colour bColour = bProgenitor->progenitor()->dataPtr()->iColour();
PDT::Colour cColour = cProgenitor->progenitor()->dataPtr()->iColour();
PDT::Colour aColour = aProgenitor->progenitor()->dataPtr()->iColour();
// decaying colour singlet
if (bColour==PDT::Colour0 ) {
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3) ||
(cColour==PDT::Colour8 && aColour==PDT::Colour8)){
dipoles.push_back(FFa);
dipoles.push_back(FFc);
}
}
// decaying colour triplet
else if (bColour==PDT::Colour3 ) {
if (cColour==PDT::Colour3 && aColour==PDT::Colour0){
dipoles.push_back(IFbc);
dipoles.push_back(IFc );
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour3){
dipoles.push_back(IFba);
dipoles.push_back(IFa );
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3){
dipoles.push_back(IFbc);
dipoles.push_back(IFc );
dipoles.push_back(FFc );
dipoles.push_back(FFa );
}
else if (cColour==PDT::Colour3 && aColour==PDT::Colour8){
dipoles.push_back(IFba);
dipoles.push_back(IFa );
dipoles.push_back(FFc );
dipoles.push_back(FFa );
}
}
// decaying colour anti-triplet
else if (bColour==PDT::Colour3bar) {
if ((cColour==PDT::Colour3bar && aColour==PDT::Colour0)){
dipoles.push_back(IFbc);
dipoles.push_back(IFc );
}
else if ((cColour==PDT::Colour0 && aColour==PDT::Colour3bar)){
dipoles.push_back(IFba);
dipoles.push_back(IFa );
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour3bar){
dipoles.push_back(IFbc);
dipoles.push_back(IFc );
dipoles.push_back(FFc );
dipoles.push_back(FFa );
}
else if (cColour==PDT::Colour3bar && aColour==PDT::Colour8){
dipoles.push_back(IFba);
dipoles.push_back(IFa );
dipoles.push_back(FFc );
dipoles.push_back(FFa );
}
}
// decaying colour octet
else if (bColour==PDT::Colour8){
if ((cColour==PDT::Colour3 && aColour==PDT::Colour3bar) ||
(cColour==PDT::Colour3bar && aColour==PDT::Colour3)){
dipoles.push_back(IFba);
dipoles.push_back(IFbc);
dipoles.push_back(IFa);
dipoles.push_back(IFc);
}
else if (cColour==PDT::Colour8 && aColour==PDT::Colour0){
dipoles.push_back(IFbc);
dipoles.push_back(IFc);
}
else if (cColour==PDT::Colour0 && aColour==PDT::Colour8){
dipoles.push_back(IFba);
dipoles.push_back(IFa);
}
}
// check colour structure is allowed
return !dipoles.empty();
}
const GeneralTwoBodyDecayer::CFlow &
GeneralTwoBodyDecayer::colourFlows(const Particle & inpart,
const ParticleVector & decay) {
// static initialization of commonly used colour structures
static const CFlow init = CFlow(3, CFlowPairVec(1, make_pair(0, 1.)));
static CFlow tripflow = init;
static CFlow atripflow = init;
static CFlow octflow = init;
static const CFlow fpflow = CFlow(4, CFlowPairVec(1, make_pair(0, 1.)));
static bool initialized = false;
if (! initialized) {
tripflow[2].resize(2, make_pair(0,1.));
tripflow[2][0] = make_pair(0, 1.);
tripflow[2][1] = make_pair(1,-1.);
tripflow[1][0] = make_pair(1, 1.);
atripflow[1].resize(2, make_pair(0,1.));
atripflow[1][0] = make_pair(0, 1.);
atripflow[1][1] = make_pair(1,-1.);
atripflow[2][0] = make_pair(1, 1.);
octflow[0].resize(2, make_pair(0,1.));
octflow[0][0] = make_pair(0,-1.);
octflow[0][1] = make_pair(1, 1.);
octflow[2][0] = make_pair(1, 1.);
initialized = true;
}
// main function body
int sing=0,trip=0,atrip=0,oct=0;
for (size_t it=0; it<decay.size(); ++it) {
switch ( decay[it]->dataPtr()->iColour() ) {
case PDT::Colour0: ++sing; break;
case PDT::Colour3: ++trip; break;
case PDT::Colour3bar: ++atrip; break;
case PDT::Colour8: ++oct; break;
}
}
// require a gluon
assert(oct>=1);
const CFlow * retval = 0;
bool inconsistent4PV = true;
// decaying colour triplet
if(inpart.dataPtr()->iColour() == PDT::Colour3 &&
trip==1 && oct==2) {
retval = &tripflow;
}
// decaying colour anti-triplet
else if(inpart.dataPtr()->iColour() == PDT::Colour3bar &&
atrip==1 && oct==2){
retval = &atripflow;
}
// decaying colour octet
else if(inpart.dataPtr()->iColour() == PDT::Colour8 &&
oct==1 && trip==1 && atrip==1) {
retval = &octflow;
}
else {
inconsistent4PV = false;
retval = &init;
}
// if a 4 point vertex exists, add a colour flow for it
if ( fourPointVertex_ ) {
if ( inconsistent4PV )
throw Exception() << "Unknown colour flows for 4 point vertex in "
<< "GeneralTwoBodyDecayer::colourFlows()"
<< Exception::runerror;
else {
retval = &fpflow;
}
}
return *retval;
}
void GeneralTwoBodyDecayer::getColourLines(vector<ColinePtr> & newline,
const HardTreePtr & hardtree,
const ShowerProgenitorPtr & bProgenitor){
// set up the colour lines
vector<ShowerParticlePtr> branchingPart;
for(set<HardBranchingPtr>::const_iterator cit=hardtree->branchings().begin();
cit!=hardtree->branchings().end();++cit)
branchingPart.push_back((**cit).branchingParticle());
static vector<int> sing,trip,atrip,oct;
sing.clear(); trip.clear(); atrip.clear(); oct.clear();
for (size_t ib=0;ib<branchingPart.size();++ib){
if (branchingPart[ib]->dataPtr()->iColour()==PDT::Colour0 ) sing. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3 ) trip. push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour3bar) atrip.push_back(ib);
else if(branchingPart[ib]->dataPtr()->iColour()==PDT::Colour8 ) oct. push_back(ib);
}
// decaying colour singlet
if (bProgenitor->progenitor()->dataPtr()->iColour()==PDT::Colour0){
if (trip.size()==1 && atrip.size()==1){
newline.push_back(new_ptr(ColourLine()));
newline[0]->addColoured (branchingPart[trip [0]]);
newline[0]->addAntiColoured(branchingPart[atrip[0]]);
}
else if (oct.size()==2){
newline.push_back(new_ptr(ColourLine()));
newline.push_back(new_ptr(ColourLine()));
newline[0]->addColoured (branchingPart[oct[0]]);
newline[0]->addAntiColoured(branchingPart[oct[1]]);
newline[1]->addColoured (branchingPart[oct[1]]);
newline[1]->addAntiColoured(branchingPart[oct[0]]);
}
}
// decaying colour triplet
else if (bProgenitor->progenitor()->dataPtr()->iColour()==PDT::Colour3 ){
if (trip.size()==2 && sing.size()==1){
newline.push_back(new_ptr(ColourLine()));
newline[0]->addColoured(branchingPart[trip[0]]);
newline[0]->addColoured(branchingPart[trip[1]]);
}
else if (trip.size()==2 && oct.size()==1){
newline.push_back(new_ptr(ColourLine()));
newline.push_back(new_ptr(ColourLine()));
if (branchingPart[trip[0]]->id()==bProgenitor->progenitor()->id()){
newline[0]->addColoured(branchingPart[trip[0]]);
newline[1]->addColoured(branchingPart[trip[1]]);
}
else {
newline[0]->addColoured(branchingPart[trip[1]]);
newline[1]->addColoured(branchingPart[trip[0]]);
}
newline[0]->addColoured (branchingPart[ oct[0]]);
newline[1]->addAntiColoured(branchingPart[ oct[0]]);
}
}
// decaying colour anti-triplet
else if (bProgenitor->progenitor()->dataPtr()->iColour()==PDT::Colour3bar){
if (atrip.size()==2 && sing.size()==1){
newline.push_back(new_ptr(ColourLine()));
newline[0]->addAntiColoured(branchingPart[atrip[0]]);
newline[0]->addAntiColoured(branchingPart[atrip[1]]);
}
else if (atrip.size()==2 && oct.size()==1){
newline.push_back(new_ptr(ColourLine()));
newline.push_back(new_ptr(ColourLine()));
if (branchingPart[atrip[0]]->id()==bProgenitor->progenitor()->id()){
newline[0]->addAntiColoured(branchingPart[atrip[0]]);
newline[1]->addAntiColoured(branchingPart[atrip[1]]);
}
else {
newline[0]->addAntiColoured(branchingPart[atrip[1]]);
newline[1]->addAntiColoured(branchingPart[atrip[0]]);
}
newline[0]->addAntiColoured(branchingPart[ oct[0]]);
newline[1]->addColoured (branchingPart[ oct[0]]);
}
}
// decaying colour octet
else if(bProgenitor->progenitor()->dataPtr()->iColour()==PDT::Colour8 ){
if (trip.size()==1 && atrip.size()==1) {
newline.push_back(new_ptr(ColourLine()));
newline.push_back(new_ptr(ColourLine()));
newline[0]->addColoured (branchingPart[ oct[0]]);
newline[1]->addAntiColoured(branchingPart[ oct[0]]);
newline[0]->addColoured (branchingPart[ trip[0]]);
newline[1]->addAntiColoured(branchingPart[atrip[0]]);
}
else if (sing.size()==1 && oct.size()==2){
newline.push_back(new_ptr(ColourLine()));
newline.push_back(new_ptr(ColourLine()));
newline[0]->addColoured (branchingPart[oct[0]]);
newline[0]->addColoured (branchingPart[oct[1]]);
newline[1]->addAntiColoured(branchingPart[oct[0]]);
newline[1]->addAntiColoured(branchingPart[oct[1]]);
}
}
}
diff --git a/Models/Feynrules/python/slha2herwig b/Models/Feynrules/python/slha2herwig
--- a/Models/Feynrules/python/slha2herwig
+++ b/Models/Feynrules/python/slha2herwig
@@ -1,93 +1,111 @@
#! /usr/bin/env python
from __future__ import division
import os, sys, argparse, re, string
comment_pat = re.compile('\s*#.*$')
block_pat = re.compile('^\s*block\s+(\w+)',flags=re.I)
decay_pat = re.compile('^\s*decay\s+(\w+)\s+([-+\.\w]+)',flags=re.I)
data_pat = re.compile('^\s*((\d+\s+)+)(-?\d\S*)\s*$')
whitespace = re.compile('\s+')
PARAMS = {}
# set up the option parser for command line input
parser = argparse.ArgumentParser(
description='Modify a ThePEG model file with parameters from a matching SLHA file.'
)
parser.add_argument(
'modelfile',
metavar='ThePEG_model',
help='ThePEG model file to use as template. Must have "# SLHA #"" annotations.'
)
parser.add_argument(
'slhafile',
metavar='SLHA_file',
help='SLHA spectrum file.'
)
+parser.add_argument(
+ '-o','--output',
+ default="FRModel_slha.model",
+ help="Name for the output file"
+)
args = parser.parse_args()
-with open(args.modelfile) as f:
- template = string.Template(f.read())
-
-
with open(args.slhafile) as f:
currentblock = None
for line in f:
line = comment_pat.sub('',line.rstrip())
if not line: continue
m = block_pat.match(line)
d = decay_pat.match(line)
if m:
currentblock = m.group(1).upper()
elif d:
currentblock = None
label = 'DECAY_%s' % d.group(1)
data = float(d.group(2))
PARAMS[label] = data
elif currentblock is not None:
d = data_pat.match(line)
if d:
index = whitespace.sub('_',d.group(1).rstrip())
try:
data = float(d.group(3))
except ValueError:
continue
label = '%s_%s' % (currentblock, index)
#if label in PARAMS:
PARAMS[label] = data
-while True:
- try:
- newcontent = template.substitute(PARAMS)
- except KeyError as e:
- problemkey = e.args[0]
- name, suffix = problemkey.rsplit('_',1)
- if suffix == 'ABS':
- mass = PARAMS[name]
- try:
- mass = mass.real
- except:
- pass
- PARAMS[problemkey] = abs(mass)
- elif suffix == 'CTAU':
- hbarc = 197.3269631e-15
- width = PARAMS[name]
- ctau = (hbarc / width) if width != 0 else 0
- PARAMS[problemkey] = ctau
- elif suffix == 'WCUT':
- width = PARAMS[name]
- wcut = 10.0 * width
- PARAMS[problemkey] = wcut
- else:
- PARAMS[problemkey] = 0.0
- else:
- break
+template = open(args.modelfile)
+output = open(args.output,'w')
-
-with open('tmp2','w') as f:
- f.write(newcontent)
- f.write('\n')
+line = template.readline()
+while line:
+ split = line.split("${")
+ if (len(split) == 1 ) :
+ output.write(line)
+ else :
+ outputLine = split[0]
+ for i in range(1,len(split)) :
+ split2 = split[i].split("}")
+ key = split2[0]
+ if key in PARAMS :
+ outputLine += str(PARAMS[key]) + split2[1]
+ else :
+ name, suffix = key.rsplit('_',1)
+ if suffix == 'ABS':
+ mass = PARAMS[name]
+ try:
+ mass = mass.real
+ except:
+ pass
+ outputLine += str(abs(mass)) + split2[1]
+ elif suffix == 'CTAU':
+ hbarc = 197.3269631e-15
+ if(name in PARAMS) :
+ width = PARAMS[name]
+ ctau = (hbarc / width) if width != 0 else 0
+ else :
+ ctau = 0
+ outputLine += str(ctau) + split2[1]
+ elif suffix == 'WCUT':
+ if(name in PARAMS) :
+ width = PARAMS[name]
+ else :
+ width = 0.
+ wcut = 10.0 * width
+ outputLine += str(wcut) + split2[1]
+ elif name == 'DECAY':
+ outputLine += str(0.) + split2[1]
+ else :
+ print 'Parameter ',key,'not set in SLHA file, keeping default'
+ outputLine = '#' + outputLine + "${" + split2[0] + "}" + split2[1]
+ output.write(outputLine)
+ line = template.readline()
+output.write("\n")
+print "Output written as : ",args.output
diff --git a/Models/Feynrules/python/ufo2herwig b/Models/Feynrules/python/ufo2herwig
--- a/Models/Feynrules/python/ufo2herwig
+++ b/Models/Feynrules/python/ufo2herwig
@@ -1,605 +1,670 @@
#! /usr/bin/env python
from __future__ import division
import os, sys, pprint, argparse, re
from string import strip, Template
# add path to the ufo conversion modules
modulepath = os.path.join("@PKGLIBDIR@",'python')
sys.path.append(modulepath)
from ufo2peg import *
# set up the option parser for command line input
parser = argparse.ArgumentParser(
description='Create Herwig++ model files from Feynrules UFO input.'
)
parser.add_argument(
'ufodir',
metavar='UFO_directory',
help='the UFO model directory'
)
parser.add_argument(
'-v', '--verbose',
action="store_true",
help="print verbose output"
)
parser.add_argument(
'-n','--name',
default="FRModel",
help="set custom nametag for the model"
)
parser.add_argument(
'--ignore-skipped',
action="store_true",
help="ignore skipped vertices and produce output anyway"
)
+parser.add_argument(
+ '--split-model',
+ action="store_true",
+ help="Split the model file into pieces to improve compilation for models with many parameters"
+)
+parser.add_argument(
+ '--no-generic-loop-vertices',
+ action="store_true",
+ help="Don't include the automatically generated generic loop vertices for h->gg and h->gamma gamma"
+)
vertex_skipped = False
args = parser.parse_args()
def should_print():
return not vertex_skipped or args.ignore_skipped
modeldir = args.ufodir.rstrip('/')
modelpath, module = os.path.split(modeldir)
if modelpath:
sys.path.append(os.path.abspath(modelpath))
FR = __import__(module)
##################################################
##################################################
# get the Model name from the arguments
modelname = args.name
libname = modelname + '.so'
# define arrays and variables
#allplist = ""
parmdecls = []
parmgetters = []
parmconstr = []
doinit = []
paramstoreplace_ = []
paramstoreplace_expressions_ = []
# get external parameters for printing
parmsubs = dict( [ (p.name, p.value)
for p in FR.all_parameters
if p.nature == 'external' ] )
# evaluate python cmath
def evaluate(x):
import cmath
return eval(x,
{'cmath':cmath,
'complexconjugate':FR.function_library.complexconjugate},
parmsubs)
## get internal params into arrays
internal = ( p for p in FR.all_parameters
if p.nature == 'internal' )
#paramstoreplaceEW_ = []
#paramstoreplaceEW_expressions_ = []
# calculate internal parameters
for p in internal:
parmsubs.update( { p.name : evaluate(p.value) } )
# if 'aS' in p.value and p.name != 'aS':
# paramstoreplace_.append(p.name)
# paramstoreplace_expressions_.append(p.value)
# if 'aEWM1' in p.value and p.name != 'aEWM1':
# paramstoreplaceEW_.append(p.name)
# paramstoreplaceEW_expressions_.append(p.value)
# more arrays used for substitution in templates
paramsforstream = []
parmmodelconstr = []
# loop over parameters and fill in template stuff according to internal/external and complex/real
# WARNING: Complex external parameter input not tested!
if args.verbose:
print 'verbose mode on: printing all parameters'
print '-'*60
paramsstuff = ('name', 'expression', 'default value', 'nature')
pprint.pprint(paramsstuff)
interfacedecl_T = """\
static Parameter<{modelname}, {type}> interface{pname}
("{pname}",
"The interface for parameter {pname}",
&{modelname}::{pname}_, {value}, 0, 0,
false, false, Interface::nolimits);
"""
interfaceDecls = []
modelparameters = {}
for p in FR.all_parameters:
value = parmsubs[p.name]
if p.type == 'real':
assert( value.imag < 1.0e-16 )
value = value.real
if p.nature == 'external':
interfaceDecls.append(
interfacedecl_T.format(modelname=modelname,
pname=p.name,
value=value,
type=typemap(p.type))
)
if hasattr(p,'lhablock'):
lhalabel = '{lhablock}_{lhacode}'.format( lhablock=p.lhablock, lhacode='_'.join(map(str,p.lhacode)) )
parmmodelconstr.append('set %s:%s ${%s}' % (modelname, p.name, lhalabel))
modelparameters[lhalabel] = value
parmsubs[p.name] = lhalabel
else:
parmmodelconstr.append('set %s:%s %s' % (modelname, p.name, value))
parmsubs[p.name] = value
parmconstr.append('%s_(%s)' % (p.name, value))
else :
parmconstr.append('%s_()' % p.name)
parmsubs[p.name] = value
elif p.type == 'complex':
value = complex(value)
if p.nature == 'external':
#
# TODO: WE DO NOT HAVE COMPLEX INTERFACES IN THEPEG (yet?)
#
# interfaceDecls.append(
# interfacedecl_T.format(modelname=modelname,
# pname=p.name,
# value='Complex(%s,%s)'%(value.real,value.imag),
# type=typemap(p.type))
# )
#
# parmmodelconstr.append('set %s:%s (%s,%s)' % (modelname, p.name, value.real, value.imag))
parmconstr.append('%s_(%s,%s)' % (p.name, value.real, value.imag))
else :
parmconstr.append('%s_(%s,%s)' % (p.name, 0.,0.))
parmsubs[p.name] = value
else:
raise Exception('Unknown data type "%s".' % p.type)
parmdecls.append(' %s %s_;' % (typemap(p.type), p.name))
parmgetters.append(' %s %s() const { return %s_; }' % (typemap(p.type),p.name, p.name))
paramsforstream.append('%s_' % p.name)
expression, symbols = 'return %s_' % p.name, None
if p.nature != 'external':
expression, symbols = py2cpp(p.value)
text = add_brackets(expression, symbols)
text=text.replace('()()','()')
#parmgetters.append(' %s %s() const { return %s; }' % (typemap(p.type),p.name, text))
doinit.append(' %s_ = %s;' % (p.name, text) )
### special treatment
# if p.name == 'aS':
# expression = '0.25 * sqr(strongCoupling(q2)) / Constants::pi'
# elif p.name == 'aEWM1':
# expression = '4.0 * Constants::pi / sqr(electroMagneticCoupling(q2))'
# elif p.name == 'Gf':
# expression = 'generator()->standardModel()->fermiConstant() * GeV2'
# elif p.name == 'MZ':
# expression = 'getParticleData(ThePEG::ParticleID::Z0)->mass() / GeV'
if args.verbose:
pprint.pprint((p.name,p.value, value, p.nature))
+paramconstructor=': '
+
+for ncount in range(0,len(parmconstr)) :
+ paramconstructor += parmconstr[ncount]
+ if(ncount != len(parmconstr) -1) :
+ paramconstructor += ','
+ if(ncount%5 == 0 ) :
+ paramconstructor += "\n"
+
+paramout=""
+paramin =""
+for ncount in range(0,len(paramsforstream)) :
+ if(ncount !=0 ) :
+ paramout += "<< " + paramsforstream[ncount]
+ paramin += ">> " + paramsforstream[ncount]
+ else :
+ paramout += paramsforstream[ncount]
+ paramin += paramsforstream[ncount]
+ if(ncount%5 == 0 ) :
+ paramout += "\n"
+ paramin += "\n"
+
parmtextsubs = { 'parmgetters' : '\n'.join(parmgetters),
'parmdecls' : '\n'.join(parmdecls),
- 'parmconstr' : ': %s' % ','.join(parmconstr),
+ 'parmconstr' : paramconstructor,
'getters' : '',
'decls' : '',
'addVertex' : '',
'doinit' : '\n'.join(doinit),
- 'ostream' : '<< '.join(paramsforstream),
- 'istream' : '>> '.join(paramsforstream),
+ 'ostream' : paramout,
+ 'istream' : paramin ,
'refs' : '',
'parmextinter': ''.join(interfaceDecls),
- 'ModelName': modelname
+ 'ModelName': modelname,
+ 'calcfunctions': ''
}
##################################################
##################################################
##################################################
# get vertex template
VERTEX = getTemplate('Vertex.cc')
VERTEXCLASS = getTemplate('Vertex_class')
VERTEXHEADER = """\
#include "ThePEG/Helicity/Vertex/{spindirectory}/{lorentztag}Vertex.h"
"""
def write_vertex_file(subs):
newname = '%s_Vertices_%03d.cc' % (subs['ModelName'],subs['vertexnumber'])
writeFile( newname, VERTEX.substitute(subs) )
if args.verbose:
print 'verbose mode on: printing all vertices'
print '-'*60
labels = ('vertex', 'particles', 'Lorentz', 'C_L', 'C_R', 'norm')
pprint.pprint(labels)
### initial pass to find global sign
def global_sign():
return 1.0
for v in FR.all_vertices:
pids = sorted([ p.pdg_code for p in v.particles ])
if pids != [-11,11,22]: continue
coup = v.couplings
assert( len(coup) == 1 )
val = coup.values()[0].value
val = evaluate(val)
assert( val.real == 0 )
return 1 if val.imag > 0 else -1
vertexclasses, vertexheaders = [], set()
ONE_EACH = True
if ONE_EACH:
all_vertices = FR.all_vertices
else:
all_vertices = collapse_vertices(FR.all_vertices)
globalsign = global_sign()
for vertexnumber,vertex in enumerate(all_vertices,1):
lorentztag = unique_lorentztag(vertex)
vertex.herwig_skip_vertex = False
# remove vertices involving ghost fields
if 'U' in lorentztag:
vertex.herwig_skip_vertex = True
continue
gsnames = ['goldstone',
'goldstoneboson',
'GoldstoneBoson']
for p in vertex.particles:
def gstest(name):
try:
return getattr(p,name)
except AttributeError:
return False
if any(map(gstest, gsnames)):
vertex.herwig_skip_vertex = True
break
if vertex.herwig_skip_vertex:
continue
lfactors = {
'FFV' : '-complex(0,1)', # ok
'VVV' : 'complex(0,1)', # changed to fix ttbar
'VVVV' : 'complex(0,1)',
'VVS' : '-complex(0,1)',
'VSS' : '-complex(0,1)', # changed to minus to fix dL ->x1 W- d
'SSS' : '-complex(0,1)', # ok
'VVSS' : '-complex(0,1)', # ok
'VVT' : 'complex(0,2)',
'VVVT' : '-complex(0,2)',
'SSSS' : '-complex(0,1)', # ok
'FFS' : '-complex(0,1)', # ok
'SST' : 'complex(0,2)',
'FFT' : '-complex(0,8)',
'FFVT' : '-complex(0,4)',
}
try:
lf = lfactors[lorentztag]
except KeyError:
msg = 'Warning: Lorentz structure {tag} ( {ps} ) in {name} ' \
'is not supported.\n'.format(tag=lorentztag, name=vertex.name,
ps=' '.join(map(str,vertex.particles)))
sys.stderr.write(msg)
vertex.herwig_skip_vertex = True
vertex_skipped = True
continue
### Particle ids
if ONE_EACH:
plistarray = [ ','.join([ str(p.pdg_code) for p in vertex.particles ]) ]
else:
plistarray = [ ','.join([ str(p.pdg_code) for p in pl ])
for pl in vertex.particles_list ]
try:
L,pos = colors(vertex)
cf = colorfactor(vertex,L,pos)
except SkipThisVertex:
vertex.herwig_skip_vertex = True
vertex_skipped = True
continue
### classname
classname = 'V_%s' % vertex.name
### parse couplings
unique_qcd = CheckUnique()
unique_qed = CheckUnique()
coup_left = []
coup_right = []
coup_norm = []
if ONE_EACH:
items = vertex.couplings.iteritems()
else:
items = vertex.couplings
try:
for (color_idx,lorentz_idx),coupling in items:
qcd, qed = qcd_qed_orders(vertex, coupling)
unique_qcd( qcd )
unique_qed( qed )
L = vertex.lorentz[lorentz_idx]
prefactors = '(%s) * (%s) * (%s)' \
% (globalsign**(len(lorentztag)-2),lf,cf[color_idx])
ordering = ''
if lorentztag in ['FFS','FFV']:
left,right = parse_lorentz(L.structure)
if left:
coup_left.append('(%s) * (%s) * (%s)' % (prefactors,left,coupling.value))
if right:
coup_right.append('(%s) * (%s) * (%s)' % (prefactors,right,coupling.value))
if lorentztag == 'FFV':
ordering = ('if(p1->id()!=%s) {Complex ltemp=left(), rtemp=right(); left(-rtemp); right(-ltemp);}'
% vertex.particles[0].pdg_code)
elif 'T' in lorentztag :
all_coup, left_coup, right_coup, ordering = \
tensorCouplings(vertex,coupling,prefactors,L,lorentztag,pos)
coup_norm += all_coup
coup_left += left_coup
coup_right += right_coup
else:
if lorentztag == 'VSS':
if L.structure == 'P(1,3) - P(1,2)':
prefactors += ' * (-1)'
ordering = 'if(p2->id()!=%s){norm(-norm());}' \
% vertex.particles[1].pdg_code
elif lorentztag == 'VVVV':
if qcd==2:
ordering = 'setType(1);\nsetOrder(0,1,2,3);'
else:
ordering, factor = EWVVVVCouplings(vertex,L)
prefactors += ' * (%s)' % factor
elif lorentztag == 'VVV':
if len(pos[8]) != 3:
ordering = VVVordering(vertex)
if type(coupling) is not list:
value = coupling.value
else:
value = "("
for coup in coupling :
value += '+(%s)' % coup.value
value +=")"
coup_norm.append('(%s) * (%s)' % (prefactors,value))
#print 'Colour :',vertex.color[color_idx]
#print 'Lorentz %s:'%L.name, L.spins, L.structure
#print 'Coupling %s:'%C.name, C.value, '\nQED=%s'%qed, 'QCD=%s'%qcd
#print '---------------'
except SkipThisVertex:
vertex.herwig_skip_vertex = True
vertex_skipped = True
continue
leftcontent = ' + '.join(coup_left) if coup_left else '0'
rightcontent = ' + '.join(coup_right) if coup_right else '0'
normcontent = ' + '.join(coup_norm) if coup_norm else '1'
# #print 'Left:',leftcontent
# #print 'Right:',rightcontent
# #print 'Norm:',normcontent
# #print '---------------'
# ### do we need left/right?
symbols = set()
if 'FF' in lorentztag and lorentztag != "FFT":
# leftcalc = aStoStrongCoup(py2cpp(leftcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
# rightcalc = aStoStrongCoup(py2cpp(rightcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
leftcalc, sym = py2cpp(leftcontent)
symbols |= sym
rightcalc, sym = py2cpp(rightcontent)
symbols |= sym
left = 'left(' + leftcalc + ');'
right = 'right(' + rightcalc + ');'
else:
left = ''
right = ''
leftcalc = ''
rightcalc = ''
# normcalc = aStoStrongCoup(py2cpp(normcontent)[0], paramstoreplace_, paramstoreplace_expressions_)
normcalc, sym = py2cpp(normcontent)
symbols |= sym
# UFO is GeV by default (?)
if lorentztag in ['VVS','SSS']:
normcalc = '(%s) * GeV / UnitRemoval::E' % normcalc
elif lorentztag in ['FFT','VVT', 'SST', 'FFVT', 'VVVT' ]:
normcalc = 'Complex((%s) / GeV * UnitRemoval::E)' % normcalc
norm = 'norm(' + normcalc + ');'
# define unkown symbols from the model
symboldefs = [ def_from_model(FR,s) for s in symbols ]
couplingptrs = [',tcPDPtr']*len(lorentztag)
if lorentztag == 'VSS':
couplingptrs[1] += ' p2'
elif lorentztag == 'FFV':
couplingptrs[0] += ' p1'
elif lorentztag == 'VVV' and ordering != '':
couplingptrs[0] += ' p1'
couplingptrs[1] += ' p2'
couplingptrs[2] += ' p3'
elif lorentztag == 'VVVT' and ordering != '':
couplingptrs[0] += ' p1'
couplingptrs[1] += ' p2'
couplingptrs[2] += ' p3'
elif lorentztag == 'VVVV' and qcd != 2:
couplingptrs[0] += ' p1'
couplingptrs[1] += ' p2'
couplingptrs[2] += ' p3'
couplingptrs[3] += ' p4'
### assemble dictionary and fill template
subs = { 'lorentztag' : lorentztag, # ok
'classname' : classname, # ok
'symbolrefs' : '\n '.join(symboldefs),
'left' : left, # doesn't always exist in base
'right' : right, # doesn't always exist in base
'norm' : norm, # needs norm, too
#################### need survey which different setter methods exist in base classes
'addToPlist' : '\n'.join([ 'addToList(%s);'%s for s in plistarray]),
'parameters' : '',
'setCouplings' : '',
'qedorder' : qed,
'qcdorder' : qcd,
'couplingptrs' : ''.join(couplingptrs),
'spindirectory' : spindirectory(lorentztag),
'ModelName' : modelname,
'ordering' : ordering,
} # ok
if args.verbose:
print '-'*60
pprint.pprint(( classname, plistarray, leftcalc, rightcalc, normcalc ))
vertexclasses.append(VERTEXCLASS.substitute(subs))
vertexheaders.add(VERTEXHEADER.format(**subs))
WRAP = 25
if vertexnumber % WRAP == 0:
write_vertex_file({'classname':classname,
'vertexnumber' : vertexnumber//WRAP,
'vertexclasses' : '\n'.join(vertexclasses),
'vertexheaders' : ''.join(vertexheaders),
'ModelName' : modelname})
vertexclasses = []
vertexheaders = set()
if not should_print():
sys.stderr.write(
"""
Error: The conversion was unsuccessful, some vertices could not be
generated. If you think the missing vertices are not important
and want to go ahead anyway, use --ignore-skipped.
Herwig++ may not give correct results, though.
"""
)
sys.exit(1)
if vertexclasses:
write_vertex_file({'classname':classname,
'vertexnumber' : vertexnumber//WRAP + 1,
'vertexclasses' : '\n'.join(vertexclasses),
'vertexheaders' : ''.join(vertexheaders),
'ModelName' : modelname})
print '='*60
##################################################
##################################################
##################################################
vertexline = """\
create Herwig::{modelname}V_{vname} /Herwig/{modelname}/V_{vname}
insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_{vname}
"""
def get_vertices():
vlist = ['library %s\n' % libname]
for v in all_vertices:
if v.herwig_skip_vertex: continue
vlist.append( vertexline.format(modelname=modelname, vname=v.name) )
- vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHPP\n'.format(modelname=modelname) )
- vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHGG\n'.format(modelname=modelname) )
+ if( not args.no_generic_loop_vertices) :
+ vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHPP\n'.format(modelname=modelname) )
+ vlist.append('insert {modelname}:ExtraVertices 0 /Herwig/{modelname}/V_GenericHGG\n'.format(modelname=modelname) )
return ''.join(vlist)
plist, names = thepeg_particles(FR,parmsubs,modelname,modelparameters)
particlelist = [
"# insert HPConstructor:Outgoing 0 /Herwig/{n}/Particles/{p}".format(n=modelname,p=p)
for p in names
]
# make the first one active to have something runnable in the example .in file
particlelist[0] = particlelist[0][2:]
particlelist = '\n'.join(particlelist)
modelfilesubs = { 'plist' : plist,
'vlist' : get_vertices(),
'setcouplings': '\n'.join(parmmodelconstr),
'ModelName': modelname
}
# write the files from templates according to the above subs
if should_print():
MODEL_HWIN = getTemplate('LHC-FR.in')
+ if(not args.split_model) :
+ MODEL_CC = [getTemplate('Model.cc')]
+ else :
+ MODEL_EXTRA_CC=getTemplate('Model6.cc')
+ extra_names=[]
+ extra_calcs=[]
+ parmtextsubs['doinit']=""
+ for i in range(0,len(doinit)) :
+ if( i%20 == 0 ) :
+ function_name = "initCalc" +str(int(i/20))
+ parmtextsubs['doinit'] += function_name +"();\n"
+ parmtextsubs['calcfunctions'] += "void " + function_name + "();\n"
+ extra_names.append(function_name)
+ extra_calcs.append("")
+ extra_calcs[-1] += doinit[i] + "\n"
+ for i in range(0,len(extra_names)) :
+ ccname = '%s_extra_%s.cc' % (modelname,i)
+ writeFile( ccname, MODEL_EXTRA_CC.substitute({'ModelName' : modelname,
+ 'functionName' : extra_names[i],
+ 'functionCalcs' : extra_calcs[i] }) )
+
+
+
+
+ MODEL_CC = [getTemplate('Model1.cc'),getTemplate('Model2.cc'),getTemplate('Model3.cc'),
+ getTemplate('Model4.cc'),getTemplate('Model5.cc')]
MODEL_H = getTemplate('Model.h')
- MODEL_CC = getTemplate('Model.cc')
+ print 'LENGTH',len(MODEL_CC)
MODELINFILE = getTemplate('FR.model')
writeFile( 'LHC-%s.in' % modelname,
MODEL_HWIN.substitute({ 'ModelName' : modelname,
'Particles' : particlelist })
)
modeltemplate = Template( MODELINFILE.substitute(modelfilesubs) )
writeFile( '%s.h' % modelname, MODEL_H.substitute(parmtextsubs) )
- writeFile( '%s.cc' % modelname, MODEL_CC.substitute(parmtextsubs) )
+ for i in range(0,len(MODEL_CC)) :
+ if(len(MODEL_CC)==1) :
+ ccname = '%s.cc' % modelname
+ else :
+ ccname = '%s.cc' % (modelname + str(i))
+ writeFile( ccname, MODEL_CC[i].substitute(parmtextsubs) )
writeFile( modelname +'.template', modeltemplate.template )
writeFile( modelname +'.model', modeltemplate.substitute( modelparameters ) )
# copy the Makefile-FR to current directory,
# replace with the modelname for compilation
with open(os.path.join(modulepath,'Makefile-FR'),'r') as orig:
with open('Makefile','w') as dest:
dest.write(orig.read().replace("FeynrulesModel.so", libname))
print 'finished generating model:\t', modelname
print 'model directory:\t\t', args.ufodir
print 'generated:\t\t\t', len(FR.all_vertices), 'vertices'
print '='*60
print 'library:\t\t\t', libname
print 'input file:\t\t\t', 'LHC-' + modelname +'.in'
print 'model file:\t\t\t', modelname +'.model'
print '='*60
print """\
To complete the installation, compile by typing "make".
An example input file is provided as LHC-FRModel.in,
you'll need to change the required particles in there.
"""
print 'DONE!'
print '='*60
diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model.cc.template
--- a/Models/Feynrules/python/ufo2peg/Model.cc.template
+++ b/Models/Feynrules/python/ufo2peg/Model.cc.template
@@ -1,72 +1,72 @@
// -*- C++ -*-
//
// ${ModelName}.cc 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ${ModelName} class.
//
#include "${ModelName}.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
//#include "Herwig++/Models/General/ModelGenerator.h"
using namespace ThePEG;
using namespace Herwig;
${ModelName}::${ModelName}()
${parmconstr}
{}
IBPtr ${ModelName}::clone() const {
return new_ptr(*this);
}
IBPtr ${ModelName}::fullclone() const {
return new_ptr(*this);
}
void ${ModelName}::doinit() {
${doinit}
- StandardModel::doinit();
+ BSMModel::doinit();
${addVertex}
}
void ${ModelName}::persistentOutput(PersistentOStream & os) const {
os << ${ostream} ;
}
void ${ModelName}::persistentInput(PersistentIStream & is, int) {
is >> ${istream} ;
}
// Static variable needed for the type description system in ThePEG.
-DescribeClass<${ModelName},StandardModel>
+DescribeClass<${ModelName},BSMModel>
describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
void ${ModelName}::Init() {
${refs}
static ClassDocumentation<${ModelName}> documentation
- ("The ${ModelName} class inherits from StandardModel"
+ ("The ${ModelName} class inherits from BSMModel"
"and supplies additional couplings and access to the ${ModelName}"
"vertices for helicity amplitude calculations" );
${parmextinter}
}
diff --git a/Models/Feynrules/python/ufo2peg/Model.h.template b/Models/Feynrules/python/ufo2peg/Model.h.template
--- a/Models/Feynrules/python/ufo2peg/Model.h.template
+++ b/Models/Feynrules/python/ufo2peg/Model.h.template
@@ -1,128 +1,137 @@
// -*- C++ -*-
//
// ${ModelName}.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2013 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_${ModelName}_H
#define HERWIG_${ModelName}_H
// This is the declaration of the ${ModelName} class.
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
namespace Herwig {
using namespace ThePEG;
using ThePEG::Constants::pi;
const Complex ii = Complex(0,1);
/** \ingroup Models
*
* This is the Herwig++ ${ModelName} class which inherits from ThePEG
* FeynRules Model class and implements additional FeynRules Model couplings,
* access to vertices for helicity amplitude calculations etc.
*
- * @see StandardModel
+ * @see BSMModel
*/
-//class ${ModelName}: public StandardModel {
-class ${ModelName}: public StandardModel {
+class ${ModelName}: public BSMModel {
public:
/// Default constructor
${ModelName}();
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);
//@}
/**
* Standard Init function used to initialize the interfaces.
*/
static void Init();
protected:
virtual bool registerDefaultVertices() const { return false; }
public:
/**
* Pointers to the objects handling the vertices.
*/
//@{
${getters}
${parmgetters}
//@}
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:
/**
* Initialize this object after the setup phase before saving and
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* Private and non-existent assignment operator.
*/
${ModelName} & operator=(const ${ModelName} &);
+
+private:
+
+ /**
+ * Helper functions for doinit
+ */
+ //@{
+
+${calcfunctions}
+ //@}
private:
/**
* Pointers to the vertices for ${ModelName} Model helicity amplitude
* calculations.
*/
//@{
${decls}
${parmdecls}
//@}
};
}
namespace ThePEG {
ThePEG_DECLARE_POINTERS(Herwig::${ModelName},Hw${ModelName}Ptr);
}
#endif /* HERWIG_${ModelName}_H */
diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model1.cc.template
copy from Models/Feynrules/python/ufo2peg/Model.cc.template
copy to Models/Feynrules/python/ufo2peg/Model1.cc.template
--- a/Models/Feynrules/python/ufo2peg/Model.cc.template
+++ b/Models/Feynrules/python/ufo2peg/Model1.cc.template
@@ -1,72 +1,38 @@
// -*- C++ -*-
//
// ${ModelName}.cc 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ${ModelName} class.
//
#include "${ModelName}.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
//#include "Herwig++/Models/General/ModelGenerator.h"
using namespace ThePEG;
using namespace Herwig;
${ModelName}::${ModelName}()
${parmconstr}
{}
IBPtr ${ModelName}::clone() const {
return new_ptr(*this);
}
IBPtr ${ModelName}::fullclone() const {
return new_ptr(*this);
}
-
-void ${ModelName}::doinit() {
- ${doinit}
- StandardModel::doinit();
- ${addVertex}
-
-}
-
-void ${ModelName}::persistentOutput(PersistentOStream & os) const {
- os << ${ostream} ;
-}
-
-void ${ModelName}::persistentInput(PersistentIStream & is, int) {
- is >> ${istream} ;
-}
-
-// Static variable needed for the type description system in ThePEG.
-DescribeClass<${ModelName},StandardModel>
- describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
-
-void ${ModelName}::Init() {
- ${refs}
-
-
-static ClassDocumentation<${ModelName}> documentation
- ("The ${ModelName} class inherits from StandardModel"
- "and supplies additional couplings and access to the ${ModelName}"
- "vertices for helicity amplitude calculations" );
-
- ${parmextinter}
-
-}
-
-
diff --git a/Models/Feynrules/python/ufo2peg/Model2.cc.template b/Models/Feynrules/python/ufo2peg/Model2.cc.template
new file mode 100644
--- /dev/null
+++ b/Models/Feynrules/python/ufo2peg/Model2.cc.template
@@ -0,0 +1,33 @@
+// -*- C++ -*-
+//
+// ${ModelName}.cc 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.
+//
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the ${ModelName} class.
+//
+
+#include "${ModelName}.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Helicity/Vertex/VertexBase.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+//#include "Herwig++/Models/General/ModelGenerator.h"
+
+using namespace ThePEG;
+using namespace Herwig;
+
+void ${ModelName}::doinit() {
+ ${doinit}
+ BSMModel::doinit();
+ ${addVertex}
+
+}
diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model3.cc.template
copy from Models/Feynrules/python/ufo2peg/Model.cc.template
copy to Models/Feynrules/python/ufo2peg/Model3.cc.template
--- a/Models/Feynrules/python/ufo2peg/Model.cc.template
+++ b/Models/Feynrules/python/ufo2peg/Model3.cc.template
@@ -1,72 +1,30 @@
// -*- C++ -*-
//
// ${ModelName}.cc 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ${ModelName} class.
//
#include "${ModelName}.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
//#include "Herwig++/Models/General/ModelGenerator.h"
using namespace ThePEG;
using namespace Herwig;
-${ModelName}::${ModelName}()
-${parmconstr}
-{}
-
-IBPtr ${ModelName}::clone() const {
- return new_ptr(*this);
-}
-
-IBPtr ${ModelName}::fullclone() const {
- return new_ptr(*this);
-}
-
-void ${ModelName}::doinit() {
- ${doinit}
- StandardModel::doinit();
- ${addVertex}
-
-}
-
void ${ModelName}::persistentOutput(PersistentOStream & os) const {
os << ${ostream} ;
}
-
-void ${ModelName}::persistentInput(PersistentIStream & is, int) {
- is >> ${istream} ;
-}
-
-// Static variable needed for the type description system in ThePEG.
-DescribeClass<${ModelName},StandardModel>
- describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
-
-void ${ModelName}::Init() {
- ${refs}
-
-
-static ClassDocumentation<${ModelName}> documentation
- ("The ${ModelName} class inherits from StandardModel"
- "and supplies additional couplings and access to the ${ModelName}"
- "vertices for helicity amplitude calculations" );
-
- ${parmextinter}
-
-}
-
-
diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model4.cc.template
copy from Models/Feynrules/python/ufo2peg/Model.cc.template
copy to Models/Feynrules/python/ufo2peg/Model4.cc.template
--- a/Models/Feynrules/python/ufo2peg/Model.cc.template
+++ b/Models/Feynrules/python/ufo2peg/Model4.cc.template
@@ -1,72 +1,30 @@
// -*- C++ -*-
//
// ${ModelName}.cc 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ${ModelName} class.
//
#include "${ModelName}.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
//#include "Herwig++/Models/General/ModelGenerator.h"
using namespace ThePEG;
using namespace Herwig;
-${ModelName}::${ModelName}()
-${parmconstr}
-{}
-
-IBPtr ${ModelName}::clone() const {
- return new_ptr(*this);
-}
-
-IBPtr ${ModelName}::fullclone() const {
- return new_ptr(*this);
-}
-
-void ${ModelName}::doinit() {
- ${doinit}
- StandardModel::doinit();
- ${addVertex}
-
-}
-
-void ${ModelName}::persistentOutput(PersistentOStream & os) const {
- os << ${ostream} ;
-}
-
void ${ModelName}::persistentInput(PersistentIStream & is, int) {
is >> ${istream} ;
}
-
-// Static variable needed for the type description system in ThePEG.
-DescribeClass<${ModelName},StandardModel>
- describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
-
-void ${ModelName}::Init() {
- ${refs}
-
-
-static ClassDocumentation<${ModelName}> documentation
- ("The ${ModelName} class inherits from StandardModel"
- "and supplies additional couplings and access to the ${ModelName}"
- "vertices for helicity amplitude calculations" );
-
- ${parmextinter}
-
-}
-
-
diff --git a/Models/Feynrules/python/ufo2peg/Model5.cc.template b/Models/Feynrules/python/ufo2peg/Model5.cc.template
new file mode 100644
--- /dev/null
+++ b/Models/Feynrules/python/ufo2peg/Model5.cc.template
@@ -0,0 +1,43 @@
+// -*- C++ -*-
+//
+// ${ModelName}.cc 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.
+//
+//
+// This is the implementation of the non-inlined, non-templated member
+// functions of the ${ModelName} class.
+//
+
+#include "${ModelName}.h"
+#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Interface/Reference.h"
+#include "ThePEG/Interface/Parameter.h"
+#include "ThePEG/Persistency/PersistentOStream.h"
+#include "ThePEG/Persistency/PersistentIStream.h"
+#include "ThePEG/Helicity/Vertex/VertexBase.h"
+#include "ThePEG/Utilities/DescribeClass.h"
+
+//#include "Herwig++/Models/General/ModelGenerator.h"
+
+using namespace ThePEG;
+using namespace Herwig;
+
+// Static variable needed for the type description system in ThePEG.
+DescribeClass<${ModelName},BSMModel>
+ describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
+
+void ${ModelName}::Init() {
+ ${refs}
+
+
+static ClassDocumentation<${ModelName}> documentation
+ ("The ${ModelName} class inherits from BSMModel"
+ "and supplies additional couplings and access to the ${ModelName}"
+ "vertices for helicity amplitude calculations" );
+
+ ${parmextinter}
+
+}
diff --git a/Models/Feynrules/python/ufo2peg/Model.cc.template b/Models/Feynrules/python/ufo2peg/Model6.cc.template
copy from Models/Feynrules/python/ufo2peg/Model.cc.template
copy to Models/Feynrules/python/ufo2peg/Model6.cc.template
--- a/Models/Feynrules/python/ufo2peg/Model.cc.template
+++ b/Models/Feynrules/python/ufo2peg/Model6.cc.template
@@ -1,72 +1,31 @@
// -*- C++ -*-
//
// ${ModelName}.cc 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ${ModelName} class.
//
#include "${ModelName}.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Utilities/DescribeClass.h"
//#include "Herwig++/Models/General/ModelGenerator.h"
using namespace ThePEG;
using namespace Herwig;
-${ModelName}::${ModelName}()
-${parmconstr}
-{}
-
-IBPtr ${ModelName}::clone() const {
- return new_ptr(*this);
+void ${ModelName}::${functionName}() {
+${functionCalcs}
}
-IBPtr ${ModelName}::fullclone() const {
- return new_ptr(*this);
-}
-
-void ${ModelName}::doinit() {
- ${doinit}
- StandardModel::doinit();
- ${addVertex}
-
-}
-
-void ${ModelName}::persistentOutput(PersistentOStream & os) const {
- os << ${ostream} ;
-}
-
-void ${ModelName}::persistentInput(PersistentIStream & is, int) {
- is >> ${istream} ;
-}
-
-// Static variable needed for the type description system in ThePEG.
-DescribeClass<${ModelName},StandardModel>
- describeThePEG${ModelName}("Herwig::${ModelName}", "${ModelName}.so");
-
-void ${ModelName}::Init() {
- ${refs}
-
-
-static ClassDocumentation<${ModelName}> documentation
- ("The ${ModelName} class inherits from StandardModel"
- "and supplies additional couplings and access to the ${ModelName}"
- "vertices for helicity amplitude calculations" );
-
- ${parmextinter}
-
-}
-
-
diff --git a/Models/General/ModelGenerator.cc b/Models/General/ModelGenerator.cc
--- a/Models/General/ModelGenerator.cc
+++ b/Models/General/ModelGenerator.cc
@@ -1,523 +1,525 @@
// -*- C++ -*-
//
// ModelGenerator.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ModelGenerator class.
//
#include "ModelGenerator.h"
-#include "ThePEG/Interface/ClassDocumentation.h"
+#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Repository/CurrentGenerator.h"
#include "BSMWidthGenerator.h"
#include "Herwig++/PDT/GenericMassGenerator.h"
#include "Herwig++/Decay/DecayIntegrator.h"
#include "ThePEG/Repository/BaseRepository.h"
using namespace Herwig;
IBPtr ModelGenerator::clone() const {
return new_ptr(*this);
}
IBPtr ModelGenerator::fullclone() const {
return new_ptr(*this);
}
void ModelGenerator::persistentOutput(PersistentOStream & os) const {
os << hardProcessConstructors_ << _theDecayConstructor << particles_
<< offshell_ << Offsel_ << BRnorm_ << twoBodyOnly_ << howOffShell_
<< Npoints_ << Iorder_ << BWshape_ << brMin_ << decayOutput_;
}
void ModelGenerator::persistentInput(PersistentIStream & is, int) {
is >> hardProcessConstructors_ >> _theDecayConstructor >> particles_
>> offshell_ >> Offsel_ >> BRnorm_ >> twoBodyOnly_ >> howOffShell_
>> Npoints_ >> Iorder_ >> BWshape_ >> brMin_ >> decayOutput_;
}
-ClassDescription<ModelGenerator> ModelGenerator::initModelGenerator;
-// Definition of the static class description member.
+// Static variable needed for the type description system in ThePEG.
+DescribeClass<ModelGenerator,Interfaced>
+describeThePEGModelGenerator("Herwig::ModelGenerator", "Herwig.so");
+
void ModelGenerator::Init() {
static ClassDocumentation<ModelGenerator> documentation
("This class controls the the use of BSM physics.",
"BSM physics was produced using the algorithm of "
"\\cite{Gigg:2007cr,Gigg:2008yc}",
"\\bibitem{Gigg:2007cr} M.~Gigg and P.~Richardson, \n"
"Eur.\\ Phys.\\ J.\\ C {\\bf 51} (2007) 989.\n"
"%%CITATION = EPHJA,C51,989;%%\n"
" %\\cite{Gigg:2008yc}\n"
"\\bibitem{Gigg:2008yc}\n"
" M.~A.~Gigg and P.~Richardson,\n"
" %``Simulation of Finite Width Effects in Physics Beyond the Standard Model,''\n"
" arXiv:0805.3037 [hep-ph].\n"
" %%CITATION = ARXIV:0805.3037;%%\n"
);
static RefVector<ModelGenerator,HardProcessConstructor>
interfaceHardProcessConstructors
("HardProcessConstructors",
"The objects to construct hard processes",
&ModelGenerator::hardProcessConstructors_, -1,
false, false, true, false, false);
static Reference<ModelGenerator,Herwig::DecayConstructor>
interfaceDecayConstructor
("DecayConstructor",
"Pointer to DecayConstructor helper class",
&ModelGenerator::_theDecayConstructor, false, false, true, false);
static RefVector<ModelGenerator,ThePEG::ParticleData> interfaceModelParticles
("DecayParticles",
"ParticleData pointers to the particles requiring spin correlation "
"decayers. If decay modes do not exist they will also be created.",
&ModelGenerator::particles_, -1, false, false, true, false);
static RefVector<ModelGenerator,ParticleData> interfaceOffshell
("Offshell",
"The particles to treat as off-shell",
&ModelGenerator::offshell_, -1, false, false, true, false);
static Switch<ModelGenerator,int> interfaceWhichOffshell
("WhichOffshell",
"A switch to determine which particles to create mass and width "
"generators for.",
&ModelGenerator::Offsel_, 0, false, false);
static SwitchOption interfaceWhichOffshellSelected
(interfaceWhichOffshell,
"Selected",
"Only create mass and width generators for the particles specified",
0);
static SwitchOption interfaceWhichOffshellAll
(interfaceWhichOffshell,
"All",
"Treat all particles specified in the DecayParticles "
"list as off-shell",
1);
static Switch<ModelGenerator,bool> interfaceBRNormalize
("BRNormalize",
"Whether to normalize the partial widths to BR*total width for an "
"on-shell particle",
&ModelGenerator::BRnorm_, true, false, false);
static SwitchOption interfaceBRNormalizeNormalize
(interfaceBRNormalize,
"Yes",
"Normalize the partial widths",
true);
static SwitchOption interfaceBRNormalizeNoNormalize
(interfaceBRNormalize,
"No",
"Do not normalize the partial widths",
false);
static Parameter<ModelGenerator,int> interfacePoints
("InterpolationPoints",
"Number of points to use for interpolation tables when needed",
&ModelGenerator::Npoints_, 10, 5, 1000,
false, false, true);
static Parameter<ModelGenerator,unsigned int>
interfaceInterpolationOrder
("InterpolationOrder", "The interpolation order for the tables",
&ModelGenerator::Iorder_, 1, 1, 5,
false, false, Interface::limited);
static Switch<ModelGenerator,int> interfaceBreitWignerShape
("BreitWignerShape",
"Controls the shape of the mass distribution generated",
&ModelGenerator::BWshape_, 0, false, false);
static SwitchOption interfaceBreitWignerShapeDefault
(interfaceBreitWignerShape,
"Default",
"Running width with q in numerator and denominator width factor",
0);
static SwitchOption interfaceBreitWignerShapeFixedWidth
(interfaceBreitWignerShape,
"FixedWidth",
"Use a fixed width",
1);
static SwitchOption interfaceBreitWignerShapeNoq
(interfaceBreitWignerShape,
"Noq",
"Use M rather than q in the numerator and denominator width factor",
2);
static SwitchOption interfaceBreitWignerShapeNoNumerator
(interfaceBreitWignerShape,
"NoNumerator",
"Neglect the numerator factors",
3);
static Switch<ModelGenerator,bool> interfaceTwoBodyOnly
("TwoBodyOnly",
"Whether to use only two-body or all modes in the running width calculation",
&ModelGenerator::twoBodyOnly_, false, false, false);
static SwitchOption interfaceTwoBodyOnlyYes
(interfaceTwoBodyOnly,
"Yes",
"Only use two-body modes",
true);
static SwitchOption interfaceTwoBodyOnlyNo
(interfaceTwoBodyOnly,
"No",
"Use all modes",
false);
static Parameter<ModelGenerator,double> interfaceMinimumBR
("MinimumBR",
"The minimum branching fraction to include",
&ModelGenerator::brMin_, 1e-6, 0.0, 1.0,
false, false, Interface::limited);
static Switch<ModelGenerator,unsigned int> interfaceDecayOutput
("DecayOutput",
"Option to control the output of the decay mode information",
&ModelGenerator::decayOutput_, 1, false, false);
static SwitchOption interfaceDecayOutputNone
(interfaceDecayOutput,
"None",
"No output",
0);
static SwitchOption interfaceDecayOutputPlain
(interfaceDecayOutput,
"Plain",
"Default plain text output",
1);
static SwitchOption interfaceDecayOutputSLHA
(interfaceDecayOutput,
"SLHA",
"Output in the Susy Les Houches Accord format",
2);
static Parameter<ModelGenerator,double> interfaceMinimumWidthFraction
("MinimumWidthFraction",
"Minimum fraction of the particle's mass the width can be"
" for the off-shell treatment.",
&ModelGenerator::minWidth_, 1e-6, 1e-15, 1.,
false, false, Interface::limited);
static Parameter<ModelGenerator,double> interfaceHowMuchOffShell
("HowMuchOffShell",
"The multiple of the particle's width by which it is allowed to be off-shell",
&ModelGenerator::howOffShell_, 5., 0.0, 100.,
false, false, Interface::limited);
}
namespace {
/// Helper function for sorting by mass
inline bool massIsLess(tcPDPtr a, tcPDPtr b) {
return a->mass() < b->mass();
}
// Helper function to find minimum possible mass of a particle
inline Energy minimumMass(tcPDPtr parent) {
Energy output(Constants::MaxEnergy);
for(set<tDMPtr>::const_iterator dit = parent->decayModes().begin();
dit != parent->decayModes().end(); ++dit) {
Energy outMass(ZERO);
for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix) {
outMass += (**dit).orderedProducts()[ix]->massMin();
}
output = min(output,outMass);
}
return output;
}
}
void ModelGenerator::doinit() {
useMe();
Interfaced::doinit();
// make sure the model is initialized
Ptr<Herwig::StandardModel>::pointer model
= dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>(generator()->standardModel());
model->init();
// and the vertices
for(size_t iv = 0; iv < model->numberOfVertices(); ++iv)
model->vertex(iv)->init();
// uniq and sort DecayParticles list by mass
set<PDPtr> tmp(particles_.begin(),particles_.end());
particles_.assign(tmp.begin(),tmp.end());
sort(particles_.begin(),particles_.end(),massIsLess);
//create decayers and decaymodes (if necessary)
if( _theDecayConstructor ) {
_theDecayConstructor->init();
_theDecayConstructor->createDecayers(particles_,brMin_);
}
// write out decays with spin correlations
ostream & os = CurrentGenerator::current().misc();
ofstream ofs;
if ( decayOutput_ > 1 ) {
string filename
= CurrentGenerator::current().filename() + "-BR.spc";
ofs.open(filename.c_str());
}
if(decayOutput_!=0) {
if(decayOutput_==1) {
os << "# The decay modes listed below will have spin\n"
<< "# correlations included when they are generated.\n#\n#";
}
else {
ofs << "# Herwig++ decay tables in SUSY Les Houches accord format\n";
ofs << "Block DCINFO # Program information\n";
ofs << "1 Herwig++ # Decay Calculator\n";
ofs << "2 " << generator()->strategy()->versionstring()
<< " # Version number\n";
}
}
//create mass and width generators for the requested particles
set<PDPtr> offShell;
if( Offsel_ == 0 ) offShell = set<PDPtr>(offshell_.begin() ,offshell_.end() );
else offShell = set<PDPtr>(particles_.begin(),particles_.end());
for(PDVector::iterator pit = particles_.begin();
pit != particles_.end(); ++pit) {
tPDPtr parent = *pit;
// Check decays for ones where quarks cannot be put on constituent
// mass-shell
checkDecays(parent);
parent->reset();
parent->update();
if( parent->CC() ) parent->CC()->synchronize();
if( parent->decaySelector().empty() ) {
parent->stable(true);
parent->width(ZERO);
parent->widthCut(ZERO);
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
else {
if(parent->mass()*minWidth_>parent->width()) {
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
else {
if( offShell.find(*pit) != offShell.end() ) {
createWidthGenerator(*pit);
}
else {
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
}
if ( decayOutput_ == 2 )
writeDecayModes(ofs, parent);
else
writeDecayModes(os, parent);
}
if( parent->massGenerator() ) {
Energy minMass = minimumMass(parent);
Energy offShellNess = howOffShell_*parent->width();
if(minMass>parent->mass()-offShellNess) {
offShellNess = parent->mass()-minMass;
}
parent->widthCut(offShellNess);
parent->massGenerator()->reset();
if(decayOutput_==1)
os << "# " <<parent->PDGName() << " will be considered off-shell.\n#\n";
}
if( parent->widthGenerator() ) parent->widthGenerator()->reset();
}
//Now construct hard processes given that we know which
//objects have running widths
for(unsigned int ix=0;ix<hardProcessConstructors_.size();++ix) {
hardProcessConstructors_[ix]->init();
hardProcessConstructors_[ix]->constructDiagrams();
}
}
void ModelGenerator::checkDecays(PDPtr parent) {
if( parent->stable() ) {
if(parent->coloured())
cerr << "Warning: No decays for coloured particle " << parent->PDGName() << "\n\n"
<< "have been calcluated in BSM model.\n"
<< "This may cause problems in the hadronization phase.\n"
<< "You may have forgotten to switch on the decay mode calculation using\n"
<< " set TwoBodyDC:CreateDecayModes Yes\n"
<< " set ThreeBodyDC:CreateDecayModes Yes\n"
<< " set WeakDecayConstructor:CreateDecayModes Yes\n"
<< "or the decays of this particle are missing from your\n"
<< "input spectrum and decay file in the SLHA format.\n\n";
return;
}
DecaySet::iterator dit = parent->decayModes().begin();
DecaySet::iterator dend = parent->decayModes().end();
Energy oldwidth(parent->width()), newwidth(ZERO);
bool rescalebrat(false);
double brsum(0.);
for(; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
Energy release((**dit).parent()->mass());
tPDVector::const_iterator pit = (**dit).orderedProducts().begin();
tPDVector::const_iterator pend =(**dit).orderedProducts().end();
for( ; pit != pend; ++pit ) {
release -= (**pit).constituentMass();
}
if( (**dit).brat() < brMin_ || release < ZERO ) {
if( release < ZERO )
cerr << "Warning: The shower cannot be generated using this decay "
<< (**dit).tag() << " because it is too close to threshold.\nIt "
<< "will be switched off and the branching fractions of the "
<< "remaining modes rescaled.\n";
rescalebrat = true;
generator()->preinitInterface(*dit, "OnOff", "set", "Off");
generator()->preinitInterface(*dit, "BranchingRatio",
"set", "0.0");
DecayIntegratorPtr decayer = dynamic_ptr_cast<DecayIntegratorPtr>((**dit).decayer());
if(decayer) {
generator()->preinitInterface(decayer->fullName(), "Initialize", "set","0");
}
}
else {
brsum += (**dit).brat();
newwidth += (**dit).brat()*oldwidth;
}
}
if( ( rescalebrat || abs(brsum - 1.) > 1e-12 ) && !parent->decayModes().empty()) {
dit = parent->decayModes().begin();
dend = parent->decayModes().end();
double factor = oldwidth/newwidth;
brsum = 0.;
for( ; dit != dend; ++dit ) {
if( !(**dit).on() ) continue;
double newbrat = ((**dit).brat())*factor;
brsum += newbrat;
ostringstream brf;
brf << setprecision(13) << newbrat;
generator()->preinitInterface(*dit, "BranchingRatio",
"set", brf.str());
}
parent->width(newwidth);
if( newwidth > ZERO ) parent->cTau(hbarc/newwidth);
}
}
namespace {
struct DecayModeOrdering {
bool operator()(tcDMPtr m1, tcDMPtr m2) {
if(m1->brat()!=m2->brat()) {
return m1->brat()>m2->brat();
}
else {
if(m1->products().size()==m2->products().size()) {
ParticleMSet::const_iterator it1=m1->products().begin();
ParticleMSet::const_iterator it2=m2->products().begin();
do {
if((**it1).id()!=(**it2).id()) {
return (**it1).id()>(**it2).id();
}
++it1;
++it2;
}
while(it1!=m1->products().end()&&
it2!=m2->products().end());
assert(false);
}
else
return m1->products().size()<m2->products().size();
}
return false;
}
};
}
void ModelGenerator::writeDecayModes(ostream & os, tcPDPtr parent) const {
if(decayOutput_==0) return;
set<tcDMPtr,DecayModeOrdering> modes;
for(set<tDMPtr>::const_iterator dit = parent->decayModes().begin();
dit != parent->decayModes().end(); ++dit) {
if((**dit).on()) modes.insert((*dit));
}
if(decayOutput_==1) {
os << " Parent: " << parent->PDGName() << " Mass (GeV): "
<< parent->mass()/GeV << " Total Width (GeV): "
<< parent->width()/GeV << endl;
os << std::left << std::setw(40) << '#'
<< std::left << std::setw(20) << "Partial Width/GeV"
<< "BR\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit)
os << std::left << std::setw(40) << (**dit).tag()
<< std::left << std::setw(20) << (**dit).brat()*parent->width()/GeV
<< (**dit).brat() << '\n';
os << "#\n#";
}
else if(decayOutput_==2) {
os << "# \t PDG \t Width\n";
os << "DECAY\t" << parent->id() << "\t" << parent->width()/GeV << "\t # " << parent->PDGName() << "\n";
for(set<tcDMPtr,DecayModeOrdering>::iterator dit=modes.begin();
dit!=modes.end();++dit) {
os << "\t" << std::left << std::setw(10)
<< (**dit).brat() << "\t" << (**dit).orderedProducts().size()
<< "\t";
for(unsigned int ix=0;ix<(**dit).orderedProducts().size();++ix)
os << std::right << std::setw(10)
<< (**dit).orderedProducts()[ix]->id() ;
for(unsigned int ix=(**dit).orderedProducts().size();ix<4;++ix)
os << "\t";
os << "# " << (**dit).tag() << "\n";
}
}
}
void ModelGenerator::createWidthGenerator(tPDPtr p) {
string wn = p->fullName() + string("-WGen");
string mn = p->fullName() + string("-MGen");
GenericMassGeneratorPtr mgen = dynamic_ptr_cast<GenericMassGeneratorPtr>
(generator()->preinitCreate("Herwig::GenericMassGenerator", mn));
BSMWidthGeneratorPtr wgen = dynamic_ptr_cast<BSMWidthGeneratorPtr>
(generator()->preinitCreate("Herwig::BSMWidthGenerator", wn));
//set the particle interface
mgen->particle(p);
wgen->particle(p);
//set the generator interfaces in the ParticleData object
generator()->preinitInterface(p, "Mass_generator","set", mn);
generator()->preinitInterface(p, "Width_generator","set", wn);
//allow the branching fraction of this particle type to vary
p->variableRatio(true);
if( p->CC() ) p->CC()->variableRatio(true);
//initialize the generators
generator()->preinitInterface(mgen, "Initialize", "set", "Yes");
generator()->preinitInterface(wgen, "Initialize", "set", "Yes");
string norm = BRnorm_ ? "Yes" : "No";
generator()->preinitInterface(wgen, "BRNormalize", "set", norm);
string twob = twoBodyOnly_ ? "Yes" : "No";
generator()->preinitInterface(wgen, "TwoBodyOnly", "set", twob);
ostringstream os;
os << Npoints_;
generator()->preinitInterface(wgen, "Points", "set", os.str());
os.str("");
os << Iorder_;
generator()->preinitInterface(wgen, "InterpolationOrder", "set",
os.str());
os.str("");
os << BWshape_;
generator()->preinitInterface(mgen, "BreitWignerShape", "set",
os.str());
}
diff --git a/Models/General/ModelGenerator.h b/Models/General/ModelGenerator.h
--- a/Models/General/ModelGenerator.h
+++ b/Models/General/ModelGenerator.h
@@ -1,252 +1,225 @@
// -*- C++ -*-
//
// ModelGenerator.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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_ModelGenerator_H
#define HERWIG_ModelGenerator_H
//
// This is the declaration of the ModelGenerator class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "DecayConstructor.h"
#include "HardProcessConstructor.h"
#include "ModelGenerator.fh"
namespace Herwig {
using namespace ThePEG;
/**
* This class is designed to store the particles in some model and
* then call the appropriate function to setup the model
*
* @see \ref ModelGeneratorInterfaces "The interfaces"
* defined for ModelGenerator.
* @see Interfaced
*/
class ModelGenerator: public Interfaced {
public:
/**
* The default constructor.
*/
ModelGenerator() : particles_(0), offshell_(0),
Offsel_(0), BRnorm_(true),
Npoints_(10), Iorder_(1),
BWshape_(0), brMin_(1e-6), twoBodyOnly_(false),
decayOutput_(1), minWidth_(1e-6),
howOffShell_(5.) {}
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();
/**
* Overloaded function from Interfaced
*/
virtual bool preInitialize() const {
return true;
}
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:
/** @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 static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ModelGenerator> initModelGenerator;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ModelGenerator & operator=(const ModelGenerator &);
private:
/**
* Check the decay modes a given particle type. This checks whether
* the decay has quarks in the final state and that they can be put on
* mass-shell during the shower.
* @param parent The parent particle
*/
void checkDecays(PDPtr parent);
/**
* Write out the spectrum of masses and decay modes
*/
void writeDecayModes(ostream & ofs, tcPDPtr parent) const;
/**
* Create mass and width generators to simulate off-shell effects
* @param p A pointer to the ParticleData object to create
* the width and mass generators for.
*/
void createWidthGenerator(tPDPtr p);
private:
/**
* Pointer to the TwoToTwoProcessConstructor
*/
vector<HPConstructorPtr> hardProcessConstructors_;
/**
* Pointer to DecayConstructor
*/
DecayConstructorPtr _theDecayConstructor;
/**
* Vector of ParticleData pointer
*/
PDVector particles_;
/** @name Width and Mass Generator variables. */
//@{
/**
* The particles to create MassGenerator and WidthGenerators
*/
PDVector offshell_;
/**
* Which particles to treat as off-shell. 1 treats all particles in
* particles_ vector as off-shell, 0 allows selection via
* offshell_ vector.
*/
int Offsel_;
/**
* Whether to normalise the partial widths to BR*Total width for
* an on-shell particle
*/
bool BRnorm_;
/**
* The number of points to include in the interpolation table
*/
int Npoints_;
/**
* The order for the interpolation
*/
unsigned int Iorder_;
/**
* The shape of the Breit-Wigner used in the mass generation
*/
int BWshape_;
/**
* The minimum branching ratio to use
*/
double brMin_;
/**
* Whether to use only two-body or all modes for running width
*/
bool twoBodyOnly_;
//@}
/**
* Option for the outputs of the decays to a file
*/
unsigned int decayOutput_;
/**
* Minimum fraction of particle's mass width can be for off-shell
* treatment
*/
double minWidth_;
/**
* How much a particle is allowed to be offshell
*/
double howOffShell_;
};
}
-#include "ThePEG/Utilities/ClassTraits.h"
-
-namespace ThePEG {
-
-/** @cond TRAITSPECIALIZATIONS */
-
-/** This template specialization informs ThePEG about the
- * base classes of ModelGenerator. */
-template <>
-struct BaseClassTrait<Herwig::ModelGenerator,1> {
- /** Typedef of the first base class of ModelGenerator. */
- typedef Interfaced NthBase;
-};
-
-/** This template specialization informs ThePEG about the name of
- * the ModelGenerator class and the shared object where it is defined. */
-template <>
-struct ClassTraits<Herwig::ModelGenerator>
- : public ClassTraitsBase<Herwig::ModelGenerator> {
- /** Return a platform-independent class name */
- static string className() { return "Herwig::ModelGenerator"; }
-};
-
-/** @endcond */
-
-}
-
#endif /* HERWIG_ModelGenerator_H */
diff --git a/Models/General/NBodyDecayConstructorBase.cc b/Models/General/NBodyDecayConstructorBase.cc
--- a/Models/General/NBodyDecayConstructorBase.cc
+++ b/Models/General/NBodyDecayConstructorBase.cc
@@ -1,630 +1,676 @@
// -*- C++ -*-
//
// NBodyDecayConstructorBase.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the NBodyDecayConstructorBase class.
//
#include "NBodyDecayConstructorBase.h"
#include "ThePEG/Repository/EventGenerator.h"
#include "ThePEG/Repository/Repository.h"
#include "ThePEG/PDT/DecayMode.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "DecayConstructor.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "ThePEG/Utilities/DescribeClass.h"
using namespace Herwig;
using namespace ThePEG;
void NBodyDecayConstructorBase::persistentOutput(PersistentOStream & os ) const {
os << init_ << iteration_ << points_ << info_ << decayConstructor_
<< createModes_ << minReleaseFraction_ << maxBoson_ << maxList_
<< removeOnShell_ << excludeEffective_ << includeTopOnShell_
<< excludedVerticesVector_ << excludedVerticesSet_
- << excludedParticlesVector_ << excludedParticlesSet_;
+ << excludedParticlesVector_ << excludedParticlesSet_
+ << removeFlavourChangingVertices_ << removeSmallVertices_
+ << minVertexNorm_;
}
void NBodyDecayConstructorBase::persistentInput(PersistentIStream & is , int) {
is >> init_ >> iteration_ >> points_ >> info_ >> decayConstructor_
>> createModes_ >> minReleaseFraction_ >> maxBoson_ >> maxList_
>> removeOnShell_ >> excludeEffective_ >> includeTopOnShell_
>> excludedVerticesVector_ >> excludedVerticesSet_
- >> excludedParticlesVector_ >> excludedParticlesSet_;
+ >> excludedParticlesVector_ >> excludedParticlesSet_
+ >> removeFlavourChangingVertices_ >> removeSmallVertices_
+ >> minVertexNorm_;
}
+// Static variable needed for the type description system in ThePEG.
+DescribeAbstractClass<NBodyDecayConstructorBase,Interfaced>
+describeThePEGNBodyDecayConstructorBase("Herwig::NBodyDecayConstructorBase", "Herwig.so");
+
AbstractClassDescription<NBodyDecayConstructorBase>
NBodyDecayConstructorBase::initNBodyDecayConstructorBase;
// Definition of the static class description member.
void NBodyDecayConstructorBase::Init() {
static ClassDocumentation<NBodyDecayConstructorBase> documentation
("The NBodyDecayConstructorBase class is the base class for the automatic"
"construction of the decay modes");
static Switch<NBodyDecayConstructorBase,bool> interfaceInitializeDecayers
("InitializeDecayers",
"Initialize new decayers",
&NBodyDecayConstructorBase::init_, true, false, false);
static SwitchOption interfaceInitializeDecayersInitializeDecayersOn
(interfaceInitializeDecayers,
"Yes",
"Initialize new decayers to find max weights",
true);
static SwitchOption interfaceInitializeDecayersoff
(interfaceInitializeDecayers,
"No",
"Use supplied weights for integration",
false);
static Parameter<NBodyDecayConstructorBase,int> interfaceInitIteration
("InitIteration",
"Number of iterations to optimise integration weights",
&NBodyDecayConstructorBase::iteration_, 1, 0, 10,
false, false, true);
static Parameter<NBodyDecayConstructorBase,int> interfaceInitPoints
("InitPoints",
"Number of points to generate when optimising integration",
&NBodyDecayConstructorBase::points_, 1000, 100, 100000000,
false, false, true);
static Switch<NBodyDecayConstructorBase,bool> interfaceOutputInfo
("OutputInfo",
"Whether to output information about the decayers",
&NBodyDecayConstructorBase::info_, false, false, false);
static SwitchOption interfaceOutputInfoOff
(interfaceOutputInfo,
"No",
"Do not output information regarding the created decayers",
false);
static SwitchOption interfaceOutputInfoOn
(interfaceOutputInfo,
"Yes",
"Output information regarding the decayers",
true);
static Switch<NBodyDecayConstructorBase,bool> interfaceCreateDecayModes
("CreateDecayModes",
"Whether to create the ThePEG::DecayMode objects as well as the decayers",
&NBodyDecayConstructorBase::createModes_, true, false, false);
static SwitchOption interfaceCreateDecayModesOn
(interfaceCreateDecayModes,
"Yes",
"Create the ThePEG::DecayMode objects",
true);
static SwitchOption interfaceCreateDecayModesOff
(interfaceCreateDecayModes,
"No",
"Only create the Decayer objects",
false);
static Switch<NBodyDecayConstructorBase,unsigned int> interfaceRemoveOnShell
("RemoveOnShell",
"Remove on-shell diagrams as should be treated as a sequence of 1->2 decays",
&NBodyDecayConstructorBase::removeOnShell_, 1, false, false);
static SwitchOption interfaceRemoveOnShellYes
(interfaceRemoveOnShell,
"Yes",
"Remove the diagrams if neither the production of decay or the intermediate"
" can happen",
1);
static SwitchOption interfaceRemoveOnShellNo
(interfaceRemoveOnShell,
"No",
"Never remove the intermediate",
0);
static SwitchOption interfaceRemoveOnShellProduction
(interfaceRemoveOnShell,
"Production",
"Remove the diagram if the on-shell production of the intermediate is allowed",
2);
static RefVector<NBodyDecayConstructorBase,VertexBase> interfaceExcludedVertices
("ExcludedVertices",
"Vertices which are not included in the three-body decayers",
&NBodyDecayConstructorBase::excludedVerticesVector_,
-1, false, false, true, true, false);
static RefVector<NBodyDecayConstructorBase,ParticleData> interfaceExcludedIntermediates
("ExcludedIntermediates",
"Excluded intermediate particles",
&NBodyDecayConstructorBase::excludedParticlesVector_,
-1, false, false, true, true, false);
static Switch<NBodyDecayConstructorBase,bool> interfaceExcludeEffectiveVertices
("ExcludeEffectiveVertices",
"Exclude effectice vertices",
&NBodyDecayConstructorBase::excludeEffective_, true, false, false);
static SwitchOption interfaceExcludeEffectiveVerticesYes
(interfaceExcludeEffectiveVertices,
"Yes",
"Exclude the effective vertices",
true);
static SwitchOption interfaceExcludeEffectiveVerticesNo
(interfaceExcludeEffectiveVertices,
"No",
"Don't exclude the effective vertices",
false);
static Parameter<NBodyDecayConstructorBase,double> interfaceMinReleaseFraction
("MinReleaseFraction",
"The minimum energy release for a three-body decay, as a "
"fraction of the parent mass.",
&NBodyDecayConstructorBase::minReleaseFraction_, 1e-3, 0.0, 1.0,
false, false, Interface::limited);
static Switch<NBodyDecayConstructorBase,unsigned int> interfaceMaximumGaugeBosons
("MaximumGaugeBosons",
"Maximum number of electroweak gauge bosons"
" to be produced as decay products",
&NBodyDecayConstructorBase::maxBoson_, 1, false, false);
static SwitchOption interfaceMaximumGaugeBosonsNone
(interfaceMaximumGaugeBosons,
"None",
"Produce no W/Zs",
0);
static SwitchOption interfaceMaximumGaugeBosonsSingle
(interfaceMaximumGaugeBosons,
"Single",
"Produce at most one W/Zs",
1);
static SwitchOption interfaceMaximumGaugeBosonsDouble
(interfaceMaximumGaugeBosons,
"Double",
"Produce at most two W/Zs",
2);
static SwitchOption interfaceMaximumGaugeBosonsTriple
(interfaceMaximumGaugeBosons,
"Triple",
"Produce at most three W/Zs",
3);
static Switch<NBodyDecayConstructorBase,unsigned int> interfaceMaximumNewParticles
("MaximumNewParticles",
"Maximum number of particles from the list of "
"decaying particles to be allowed as decay products",
&NBodyDecayConstructorBase::maxList_, 0, false, false);
static SwitchOption interfaceMaximumNewParticlesNone
(interfaceMaximumNewParticles,
"None",
"No particles from the list",
0);
static SwitchOption interfaceMaximumNewParticlesOne
(interfaceMaximumNewParticles,
"One",
"A single particle from the list",
1);
static SwitchOption interfaceMaximumNewParticlesTwo
(interfaceMaximumNewParticles,
"Two",
"Two particles from the list",
2);
static SwitchOption interfaceMaximumNewParticlesThree
(interfaceMaximumNewParticles,
"Three",
"Three particles from the list",
3);
static SwitchOption interfaceMaximumNewParticlesFour
(interfaceMaximumNewParticles,
"Four",
"Four particles from the list",
4);
static Switch<NBodyDecayConstructorBase,bool> interfaceIncludeOnShellTop
("IncludeOnShellTop",
"Include the on-shell diagrams involving t -> bW",
&NBodyDecayConstructorBase::includeTopOnShell_, false, false, false);
static SwitchOption interfaceIncludeOnShellTopYes
(interfaceIncludeOnShellTop,
"Yes",
"Inlude them",
true);
static SwitchOption interfaceIncludeOnShellTopNo
(interfaceIncludeOnShellTop,
"No",
"Don't include them",
true);
+
+ static Switch<NBodyDecayConstructorBase,bool> interfaceRemoveSmallVertices
+ ("RemoveSmallVertices",
+ "Remove vertices with norm() below minVertexNorm",
+ &NBodyDecayConstructorBase::removeSmallVertices_, false, false, false);
+ static SwitchOption interfaceRemoveSmallVerticesYes
+ (interfaceRemoveSmallVertices,
+ "Yes",
+ "Remove them",
+ true);
+ static SwitchOption interfaceRemoveSmallVerticesNo
+ (interfaceRemoveSmallVertices,
+ "No",
+ "Don't remove them",
+ false);
+
+ static Parameter<NBodyDecayConstructorBase,double> interfaceMinVertexNorm
+ ("MinVertexNorm",
+ "Minimum allowed value of the notm() of the vertex if removing small vertices",
+ &NBodyDecayConstructorBase::minVertexNorm_, 1e-8, 1e-300, 1.,
+ false, false, Interface::limited);
+
+ static Switch<NBodyDecayConstructorBase,bool> interfaceRemoveFlavourChangingVertices
+ ("RemoveFlavourChangingVertices",
+ "Remove flavour changing interactions with the photon and gluon",
+ &NBodyDecayConstructorBase::removeFlavourChangingVertices_, false, false, false);
+ static SwitchOption interfaceRemoveFlavourChangingVerticesYes
+ (interfaceRemoveFlavourChangingVertices,
+ "Yes",
+ "Remove them",
+ true);
+ static SwitchOption interfaceRemoveFlavourChangingVerticesNo
+ (interfaceRemoveFlavourChangingVertices,
+ "No",
+ "Don't remove them",
+ false);
+
}
void NBodyDecayConstructorBase::setBranchingRatio(tDMPtr dm, Energy pwidth) {
//Need width and branching ratios for all currently created decay modes
PDPtr parent = const_ptr_cast<PDPtr>(dm->parent());
DecaySet modes = parent->decayModes();
if( modes.empty() ) return;
double dmbrat(0.);
if( modes.size() == 1 ) {
parent->width(pwidth);
if( pwidth > ZERO ) parent->cTau(hbarc/pwidth);
dmbrat = 1.;
}
else {
Energy currentwidth(parent->width());
Energy newWidth(currentwidth + pwidth);
parent->width(newWidth);
if( newWidth > ZERO ) parent->cTau(hbarc/newWidth);
//need to reweight current branching fractions if there are any
double factor = newWidth > ZERO ? double(currentwidth/newWidth) : 0.;
for(DecaySet::const_iterator dit = modes.begin();
dit != modes.end(); ++dit) {
if( **dit == *dm || !(**dit).on() ) continue;
double newbrat = (**dit).brat()*factor;
ostringstream brf;
brf << setprecision(13)<< newbrat;
generator()->preinitInterface(*dit, "BranchingRatio",
"set", brf.str());
}
//set brat for current mode
dmbrat = newWidth > ZERO ? double(pwidth/newWidth) : 0.;
}
ostringstream br;
br << setprecision(13) << dmbrat;
generator()->preinitInterface(dm, "BranchingRatio",
"set", br.str());
}
void NBodyDecayConstructorBase::setDecayerInterfaces(string fullname) const {
if( initialize() ) {
ostringstream value;
value << initialize();
generator()->preinitInterface(fullname, "Initialize", "set",
value.str());
value.str("");
value << iteration();
generator()->preinitInterface(fullname, "Iteration", "set",
value.str());
value.str("");
value << points();
generator()->preinitInterface(fullname, "Points", "set",
value.str());
}
// QED stuff if needed
if(decayConstructor()->QEDGenerator())
generator()->preinitInterface(fullname, "PhotonGenerator", "set",
decayConstructor()->QEDGenerator()->fullName());
string outputmodes;
if( info() ) outputmodes = string("Output");
else outputmodes = string("NoOutput");
generator()->preinitInterface(fullname, "OutputModes", "set",
outputmodes);
}
void NBodyDecayConstructorBase::doinit() {
Interfaced::doinit();
excludedVerticesSet_ = set<VertexBasePtr>(excludedVerticesVector_.begin(),
excludedVerticesVector_.end());
excludedParticlesSet_ = set<PDPtr>(excludedParticlesVector_.begin(),
excludedParticlesVector_.end());
if(removeOnShell_==0&&numBodies()>2)
generator()->log() << "Warning: Including diagrams with on-shell "
<< "intermediates in " << numBodies() << "-body BSM decays, this"
<< " can lead to double counting and is not"
<< " recommended unless you really know what you are doing\n"
<< "This can be switched off using\n set "
<< fullName() << ":RemoveOnShell Yes\n";
}
namespace {
double factorial(const int i) {
if(i>1) return i*factorial(i-1);
else return 1.;
}
void constructIdenticalSwaps(unsigned int depth,
vector<vector<unsigned int> > identical,
vector<unsigned int> channelType,
list<vector<unsigned int> > & swaps) {
if(depth==0) {
unsigned int size = identical.size();
for(unsigned ix=0;ix<size;++ix) {
for(unsigned int iy=2;iy<identical[ix].size();++iy)
identical.push_back(identical[ix]);
}
}
if(depth+1!=identical.size()) {
constructIdenticalSwaps(depth+1,identical,channelType,swaps);
}
else {
swaps.push_back(channelType);
}
for(unsigned int ix=0;ix<identical[depth].size();++ix) {
for(unsigned int iy=ix+1;iy<identical[depth].size();++iy) {
vector<unsigned int> newType=channelType;
swap(newType[identical[depth][ix]],newType[identical[depth][iy]]);
// at bottom of chain
if(depth+1==identical.size()) {
swaps.push_back(newType);
}
else {
constructIdenticalSwaps(depth+1,identical,newType,swaps);
}
}
}
}
void identicalFromSameDecay(unsigned int & loc, const NBVertex & vertex,
vector<vector<unsigned int> > & sameDecay) {
list<pair<PDPtr,NBVertex> >::const_iterator it = vertex.vertices.begin();
while(it!=vertex.vertices.end()) {
if(it->second.incoming) {
identicalFromSameDecay(loc,it->second,sameDecay);
++it;
continue;
}
++loc;
long id = it->first->id();
++it;
if(it == vertex.vertices.end()) break;
if(it->second.incoming) continue;
if(it->first->id()!=id) continue;
sameDecay.push_back(vector<unsigned int>());
sameDecay.back().push_back(loc-1);
while(it != vertex.vertices.end()
&& !it->second.incoming
&& it->first->id()==id) {
++loc;
++it;
sameDecay.back().push_back(loc-1);
}
};
}
}
void NBodyDecayConstructorBase::DecayList(const set<PDPtr> & particles) {
if( particles.empty() ) return;
// cast the StandardModel to the Hw++ one to get the vertices
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
unsigned int nv(model->numberOfVertices());
// loop over the particles and create the modes
for(set<PDPtr>::const_iterator ip=particles.begin();
ip!=particles.end();++ip) {
// get the decaying particle
tPDPtr parent = *ip;
if ( Debug::level > 0 )
Repository::cout() << "Constructing N-body decays for "
<< parent->PDGName() << '\n';
// first create prototype 1->2 decays
std::stack<PrototypeVertexPtr> prototypes;
for(unsigned int iv = 0; iv < nv; ++iv) {
VertexBasePtr vertex = model->vertex(iv);
if(excluded(vertex)) continue;
- PrototypeVertex::createPrototypes(parent, vertex, prototypes);
+ PrototypeVertex::createPrototypes(parent, vertex, prototypes,this);
}
// now expand them into potential decay modes
list<vector<PrototypeVertexPtr> > modes;
while(!prototypes.empty()) {
// get the first prototype from the stack
PrototypeVertexPtr proto = prototypes.top();
prototypes.pop();
// multiplcity too low
if(proto->npart<numBodies()) {
// loop over all vertices and expand
for(unsigned int iv = 0; iv < nv; ++iv) {
VertexBasePtr vertex = model->vertex(iv);
if(excluded(vertex) ||
proto->npart+vertex->getNpoint()>numBodies()+2) continue;
PrototypeVertex::expandPrototypes(proto,vertex,prototypes,
- excludedParticlesSet_);
+ excludedParticlesSet_,this);
}
}
// multiplcity too high disgard
else if(proto->npart>numBodies()) {
continue;
}
// right multiplicity
else {
// check it's kinematical allowed for physical masses
if( proto->incomingMass() < proto->outgoingMass() ) continue;
// and for constituent masses
Energy outgoingMass = proto->outgoingConstituentMass();
if( proto->incomingMass() < proto->outgoingConstituentMass() ) continue;
// remove processes which aren't kinematically allowed within
// the release fraction
if( proto->incomingMass() - outgoingMass <=
minReleaseFraction_ * proto->incomingMass() ) continue;
// check the external particles
if(!proto->checkExternal()) continue;
// check if first piece on-shell
bool onShell = proto->canBeOnShell(removeOnShell(),proto->incomingMass(),true);
// special treatment for three-body top decays
if(onShell) {
if(includeTopOnShell_ &&
abs(proto->incoming->id())==ParticleID::t && proto->npart==3) {
unsigned int nprimary=0;
bool foundW=false, foundQ=false;
for(OrderedVertices::const_iterator it = proto->outgoing.begin();
it!=proto->outgoing.end();++it) {
if(abs(it->first->id())==ParticleID::Wplus)
foundW = true;
if(abs(it->first->id())==ParticleID::b ||
abs(it->first->id())==ParticleID::s ||
abs(it->first->id())==ParticleID::d)
foundQ = true;
++nprimary;
}
if(!(foundW&&foundQ&&nprimary==2)) continue;
}
else continue;
}
// check if should be added to an existing decaymode
bool added = false;
for(list<vector<PrototypeVertexPtr> >::iterator it=modes.begin();
it!=modes.end();++it) {
// is the same ?
if(!(*it)[0]->sameDecay(*proto)) continue;
// it is the same
added = true;
// check if it is a duplicate
bool already = false;
for(unsigned int iz = 0; iz < (*it).size(); ++iz) {
if( *(*it)[iz] == *proto) {
already = true;
break;
}
}
if(!already) (*it).push_back(proto);
break;
}
if(!added) modes.push_back(vector<PrototypeVertexPtr>(1,proto));
}
}
// now look at the decay modes
for(list<vector<PrototypeVertexPtr> >::iterator mit = modes.begin();
mit!=modes.end();++mit) {
// count the number of gauge bosons and particles from the list
unsigned int nlist(0),nbos(0);
for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin();
it!=(*mit)[0]->outPart.end();++it) {
if(abs((**it).id()) == ParticleID::Wplus ||
abs((**it).id()) == ParticleID::Z0) ++nbos;
if(particles.find(*it)!=particles.end()) ++nlist;
if((**it).CC() && particles.find((**it).CC())!=particles.end()) ++nlist;
}
// if too many ignore the mode
if(nbos > maxBoson_ || nlist > maxList_) continue;
// translate the prototypes into diagrams
vector<NBDiagram> newDiagrams;
double symfac(1.);
bool possibleOnShell=false;
for(unsigned int ix=0;ix<(*mit).size();++ix) {
symfac = 1.;
possibleOnShell |= (*mit)[ix]->possibleOnShell;
NBDiagram templateDiagram = NBDiagram((*mit)[ix]);
// extract the ordering
vector<int> order(templateDiagram.channelType.size());
for(unsigned int iz=0;iz<order.size();++iz) {
order[templateDiagram.channelType[iz]-1]=iz;
}
// find any identical particles
vector<vector<unsigned int> > identical;
OrderedParticles::const_iterator it=templateDiagram.outgoing.begin();
unsigned int iloc=0;
while(it!=templateDiagram.outgoing.end()) {
OrderedParticles::const_iterator lt = templateDiagram.outgoing.lower_bound(*it);
OrderedParticles::const_iterator ut = templateDiagram.outgoing.upper_bound(*it);
unsigned int nx=0;
for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {++nx;}
if(nx==1) {
++it;
++iloc;
}
else {
symfac *= factorial(nx);
identical.push_back(vector<unsigned int>());
for(OrderedParticles::const_iterator jt=lt;jt!=ut;++jt) {
identical.back().push_back(order[iloc]);
++iloc;
}
it = ut;
}
}
// that's it if there outgoing are unqiue
if(identical.empty()) {
newDiagrams.push_back(templateDiagram);
continue;
}
// find any identical particles which shouldn't be swapped
unsigned int loc=0;
vector<vector<unsigned int> > sameDecay;
identicalFromSameDecay(loc,templateDiagram,sameDecay);
// compute the swaps
list<vector<unsigned int> > swaps;
constructIdenticalSwaps(0,identical,templateDiagram.channelType,swaps);
// special if identical from same decay
if(!sameDecay.empty()) {
for(vector<vector<unsigned int> >::const_iterator st=sameDecay.begin();
st!=sameDecay.end();++st) {
for(list<vector<unsigned int> >::iterator it=swaps.begin();
it!=swaps.end();++it) {
for(unsigned int iy=0;iy<(*st).size();++iy) {
for(unsigned int iz=iy+1;iz<(*st).size();++iz) {
if((*it)[(*st)[iy]]>(*it)[(*st)[iz]])
swap((*it)[(*st)[iy]],(*it)[(*st)[iz]]);
}
}
}
}
}
// remove any dupliciates
for(list<vector<unsigned int> >::iterator it=swaps.begin();
it!=swaps.end();++it) {
for(list<vector<unsigned int> >::iterator jt=it;
jt!=swaps.end();++jt) {
if(it==jt) continue;
bool different=false;
for(unsigned int iz=0;iz<(*it).size();++iz) {
if((*it)[iz]!=(*jt)[iz]) {
different = true;
break;
}
}
if(!different) {
jt = swaps.erase(jt);
--jt;
}
}
}
// special for identical decay products
if(templateDiagram.vertices.begin()->second.incoming) {
if( templateDiagram.vertices.begin() ->first ==
(++templateDiagram.vertices.begin())->first) {
if(*( (*mit)[ix]->outgoing.begin() ->second) ==
*((++(*mit)[ix]->outgoing.begin())->second)) {
for(list<vector<unsigned int> >::iterator it=swaps.begin();
it!=swaps.end();++it) {
for(list<vector<unsigned int> >::iterator jt=it;
jt!=swaps.end();++jt) {
if(it==jt) continue;
if((*it)[0]==(*jt)[2]&&(*it)[1]==(*jt)[3]) {
jt = swaps.erase(jt);
--jt;
}
}
}
}
}
}
for(list<vector<unsigned int> >::iterator it=swaps.begin();
it!=swaps.end();++it) {
newDiagrams.push_back(templateDiagram);
newDiagrams.back().channelType = *it;
}
}
// create the decay
if( Debug::level > 1 ) {
generator()->log() << "Mode: ";
generator()->log() << (*mit)[0]->incoming->PDGName() << " -> ";
for(OrderedParticles::const_iterator it=(*mit)[0]->outPart.begin();
it!=(*mit)[0]->outPart.end();++it)
generator()->log() << (**it).PDGName() << " ";
generator()->log() << "\n";
generator()->log() << "There are " << (*mit).size() << " diagrams\n";
for(unsigned int iy=0;iy<newDiagrams.size();++iy) {
generator()->log() << "Diagram: " << iy << "\n";
generator()->log() << newDiagrams[iy] << "\n";
}
}
createDecayMode(newDiagrams,possibleOnShell,symfac);
}
}
}
void NBodyDecayConstructorBase::createDecayMode(vector<NBDiagram> &,
bool, double) {
throw Exception() << "In NBodyDecayConstructorBase::createDecayMode() which"
<< " should have be overridden in an inheriting class"
<< Exception::abortnow;
}
diff --git a/Models/General/NBodyDecayConstructorBase.h b/Models/General/NBodyDecayConstructorBase.h
--- a/Models/General/NBodyDecayConstructorBase.h
+++ b/Models/General/NBodyDecayConstructorBase.h
@@ -1,336 +1,347 @@
// -*- C++ -*-
//
// NBodyDecayConstructorBase.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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_NBodyDecayConstructorBase_H
#define HERWIG_NBodyDecayConstructorBase_H
//
// This is the declaration of the NBodyDecayConstructorBase class.
//
#include "ThePEG/Interface/Interfaced.h"
#include "ThePEG/Utilities/Exception.h"
#include "ThePEG/PDT/ParticleData.h"
#include "NBodyDecayConstructorBase.fh"
#include "PrototypeVertex.h"
#include "DecayConstructor.fh"
namespace Herwig {
using namespace ThePEG;
/**
* This is the base class for NBodyDecayConstructors. An N-body
* decay constructor should inherit from this and implement the
* DecayList virtual funtcion to create the decays and decayers.
*
* @see \ref NBodyDecayConstructorBaseInterfaces "The interfaces"
* defined for NBodyDecayConstructor.
*/
class NBodyDecayConstructorBase: public Interfaced {
public:
/**
* The default constructor.
*/
NBodyDecayConstructorBase() :
init_(true),iteration_(1), points_(1000), info_(false),
createModes_(true), removeOnShell_(1), excludeEffective_(true),
minReleaseFraction_(1e-3), maxBoson_(1), maxList_(1),
- includeTopOnShell_(false) {}
+ includeTopOnShell_(false ), removeFlavourChangingVertices_(false),
+ removeSmallVertices_(false), minVertexNorm_(1e-8)
+ {}
/**
* 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 = 0;
/**
* Set the pointer to the DecayConstrcutor
*/
void decayConstructor(tDecayConstructorPtr d) {
decayConstructor_ = d;
}
+ /**
+ * Remove flavour changing vertices ?
+ */
+ bool removeFlavourChangingVertices() const {
+ return removeFlavourChangingVertices_;
+ }
+
+ /**
+ * Remove small vertices ?
+ */
+ bool removeSmallVertices() const {
+ return removeSmallVertices_;
+ }
+
+ /**
+ * Minimum norm for vertex removal
+ */
+ double minVertexNorm() const {
+ return minVertexNorm_;
+ }
+
protected:
/**
* Method to set up the decay mode, should be overidden in inheriting class
*/
virtual void createDecayMode(vector<NBDiagram> & mode,
bool possibleOnShell,
double symfac);
/**
* Set the branching ratio of this mode. This requires
* calculating a new width for the decaying particle and reweighting
* the current branching fractions.
* @param dm The decaymode for which to set the branching ratio
* @param pwidth The calculated width of the mode
*/
void setBranchingRatio(tDMPtr dm, Energy pwidth);
/**
* Set the interfaces of the decayers depending on the flags stored.
* @param name Fullname of the decayer in the EventGenerator
* including the path
*/
void setDecayerInterfaces(string name) const;
/**
* Whether to initialize decayers or not
*/
bool initialize() const { return init_; }
/**
* Number of iterations if initializing (default 1)
*/
int iteration() const { return iteration_; }
/**
* Number of points to do in initialization
*/
int points() const { return points_; }
/**
* Whether to output information on the decayers
*/
bool info() const { return info_; }
/**
* Whether to create the DecayModes as well as the Decayer objects
*/
bool createDecayModes() const { return createModes_; }
/**
* Maximum number of electroweak gauge bosons
*/
unsigned int maximumGaugeBosons() const { return maxBoson_;}
/**
* Maximum number of particles from the list whose decays we are calculating
*/
unsigned int maximumList() const { return maxList_;}
/**
* Minimum energy release fraction
*/
double minimumReleaseFraction() const {return minReleaseFraction_;}
/**
* Get the pointer to the DecayConstructor object
*/
tDecayConstructorPtr decayConstructor() const {
return decayConstructor_;
}
/**
* Option for on-shell particles
*/
unsigned int removeOnShell() const { return removeOnShell_; }
/**
* Check if a vertex is excluded
*/
bool excluded(VertexBasePtr vertex) const {
// skip an effective vertex
if( excludeEffective_ &&
int(vertex->orderInGs() + vertex->orderInGem()) != int(vertex->getNpoint())-2)
return true;
// check if explicitly forbidden
return excludedVerticesSet_.find(vertex)!=excludedVerticesSet_.end();
}
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();
//@}
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is an abstract class with persistent data.
*/
static AbstractClassDescription<NBodyDecayConstructorBase>
initNBodyDecayConstructorBase;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
NBodyDecayConstructorBase & operator=(const NBodyDecayConstructorBase &);
private:
/**
* Whether to initialize decayers or not
*/
bool init_;
/**
* Number of iterations if initializing (default 1)
*/
int iteration_;
/**
* Number of points to do in initialization
*/
int points_;
/**
* Whether to output information on the decayers
*/
bool info_;
/**
* Whether to create the DecayModes as well as the Decayer objects
*/
bool createModes_;
/**
* Whether or not to remove on-shell diagrams
*/
unsigned int removeOnShell_;
/**
* Excluded Vertices
*/
vector<VertexBasePtr> excludedVerticesVector_;
/**
* Excluded Vertices
*/
set<VertexBasePtr> excludedVerticesSet_;
/**
* Excluded Particles
*/
vector<PDPtr> excludedParticlesVector_;
/**
* Excluded Particles
*/
set<PDPtr> excludedParticlesSet_;
/**
* Whether or not to exclude effective vertices
*/
bool excludeEffective_;
/**
* A pointer to the DecayConstructor object
*/
tDecayConstructorPtr decayConstructor_;
/**
* The minimum energy release for a three-body decay as a
* fraction of the parent mass
*/
double minReleaseFraction_;
/**
* Maximum number of EW gauge bosons
*/
unsigned int maxBoson_;
/**
* Maximum number of particles from the decaying particle list
*/
unsigned int maxList_;
/**
* Include on-shell for \f$t\to b W\f$
*/
bool includeTopOnShell_;
+
+ /**
+ * Remove flavour changing vertices ?
+ */
+ bool removeFlavourChangingVertices_;
+
+ /**
+ * Remove small vertices ?
+ */
+ bool removeSmallVertices_;
+
+ /**
+ * Minimum norm for vertex removal
+ */
+ double minVertexNorm_;
};
/** An Exception class that can be used by all inheriting classes to
* indicate a setup problem. */
class NBodyDecayConstructorError : public Exception {
public:
NBodyDecayConstructorError() : Exception() {}
NBodyDecayConstructorError(const string & str,
Severity sev) : Exception(str,sev)
{}
};
}
-#include "ThePEG/Utilities/ClassTraits.h"
-
-namespace ThePEG {
-
-/** @cond TRAITSPECIALIZATIONS */
-
-/** This template specialization informs ThePEG about the
- * base classes of NBodyDecayConstructorBase. */
-template <>
-struct BaseClassTrait<Herwig::NBodyDecayConstructorBase,1> {
- /** Typedef of the first base class of NBodyDecayConstructorBase. */
- typedef Interfaced NthBase;
-};
-
-/** This template specialization informs ThePEG about the name of
- * the NBodyDecayConstructorBase class and the shared object where it is defined. */
-template <>
-struct ClassTraits<Herwig::NBodyDecayConstructorBase>
- : public ClassTraitsBase<Herwig::NBodyDecayConstructorBase> {
- /** Return a platform-independent class name */
- static string className() { return "Herwig::NBodyDecayConstructorBase"; }
-};
-
-/** @endcond */
-
-}
-
#endif /* HERWIG_NBodyDecayConstructorBase_H */
diff --git a/Models/General/PrototypeVertex.cc b/Models/General/PrototypeVertex.cc
--- a/Models/General/PrototypeVertex.cc
+++ b/Models/General/PrototypeVertex.cc
@@ -1,288 +1,320 @@
// -*- 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"
+#include "NBodyDecayConstructorBase.h"
using namespace Herwig;
void PrototypeVertex::createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
- std::stack<PrototypeVertexPtr> & prototypes) {
+ std::stack<PrototypeVertexPtr> & prototypes,
+ NBodyDecayConstructorBasePtr decayCon) {
if(!vertex->isIncoming(inpart)) return;
tPDPtr ccpart = inpart->CC() ? inpart->CC() : inpart;
long id = ccpart->id();
+ Energy2 scale = sqr(inpart->mass());
for(unsigned int list=0;list<vertex->getNpoint();++list) {
tPDVector decaylist = vertex->search(list, ccpart);
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]);
out.insert(make_pair(pout[ix],PrototypeVertexPtr()));
}
+ // remove flavour changing vertices
+ if(decayCon->removeFlavourChangingVertices()&&vertex->getNpoint()==3) {
+ if((pout[1]->id()==ParticleID::gamma||
+ pout[1]->id()==ParticleID::g) &&
+ inpart!=pout[2]) continue;
+ else if((pout[2]->id()==ParticleID::gamma||
+ pout[2]->id()==ParticleID::g) &&
+ inpart!=pout[1]) continue;
+ }
+ // remove small vertices if needed
+ if(decayCon->removeSmallVertices()) {
+ if(vertex->getNpoint()==3)
+ vertex->setCoupling(scale,pout[0],pout[1],pout[2]);
+ else if(vertex->getNpoint()==4)
+ vertex->setCoupling(scale,pout[0],pout[1],pout[2],pout[3]);
+ if(abs(vertex->norm())<decayCon->minVertexNorm()) continue;
+ }
if(vertex->getNpoint()==3) {
// remove radiation
- if(
- (inpart==pout[1] && (pout[2]->id()==ParticleID::gamma||
+ if((inpart==pout[1] && (pout[2]->id()==ParticleID::gamma||
pout[2]->id()==ParticleID::g||
pout[2]->id()==ParticleID::Z0))
||
(inpart==pout[2] && (pout[1]->id()==ParticleID::gamma||
pout[1]->id()==ParticleID::g||
- pout[1]->id()==ParticleID::Z0))
- ||
- (inpart->id()==ParticleID::g && pout[1]->CC() && pout[1]->CC() == pout[2])
- )
+ 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::stack<PrototypeVertexPtr> & prototypes,
- const set<PDPtr> & excluded) {
+ const set<PDPtr> & excluded,
+ NBodyDecayConstructorBasePtr decayCon) {
for(OrderedVertices::const_iterator it = proto->outgoing.begin();
it!=proto->outgoing.end();++it) {
if(it->second) {
- expandPrototypes(it->second,vertex,prototypes,excluded);
+ expandPrototypes(it->second,vertex,prototypes,excluded,decayCon);
}
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;
tcPDPtr ccpart = it->first->CC() ? it->first->CC() : it->first;
long id = ccpart->id();
PrototypeVertexPtr parent=proto;
while(parent->parent) parent=parent->parent;
+ Energy2 scale(sqr(parent->incoming->mass()));
for(unsigned int il = 0; il < vertex->getNpoint(); ++il) {
tPDVector decaylist = vertex->search(il,ccpart);
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]);
outgoing.insert(make_pair(pout[iy],PrototypeVertexPtr()));
}
if(vertex->getNpoint()==3) {
if((it->first==pout[1] && (pout[2]->id()==ParticleID::gamma||
pout[2]->id()==ParticleID::g||
pout[2]->id()==ParticleID::Z0))
||
(it->first==pout[2] && (pout[1]->id()==ParticleID::gamma||
pout[1]->id()==ParticleID::g||
- pout[1]->id()==ParticleID::Z0))
- ||
- (it->first->id()==ParticleID::g && pout[1]->CC() && pout[1]->CC() == pout[2])
- )
+ 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;
+ }
+ // remove flavour changing vertices
+ if(decayCon->removeFlavourChangingVertices()&&vertex->getNpoint()==3) {
+ if((pout[1]->id()==ParticleID::gamma||
+ pout[1]->id()==ParticleID::g) &&
+ it->first!=pout[2]) continue;
+ else if((pout[2]->id()==ParticleID::gamma||
+ pout[2]->id()==ParticleID::g) &&
+ it->first!=pout[1]) continue;
+ }
+ // remove small vertices if needed
+ if(decayCon->removeSmallVertices()) {
+ if(vertex->getNpoint()==3)
+ vertex->setCoupling(scale,pout[0],pout[1],pout[2]);
+ else if(vertex->getNpoint()==4)
+ vertex->setCoupling(scale,pout[0],pout[1],pout[2],pout[3]);
+ if(abs(vertex->norm())<decayCon->minVertexNorm()) 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;
}
bool PrototypeVertex::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;
}
namespace {
void constructTag(vector<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(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.push_back(ix+1);
break;
}
++ix;
}
}
}
}
}
NBDiagram::NBDiagram(PrototypeVertexPtr proto)
: NBVertex(proto),
colourFlow (vector<CFPair>(1,make_pair(1,1.))),
largeNcColourFlow(vector<CFPair>(1,make_pair(1,1.))) {
if(!proto) return;
// finally work out the channel and the ordering
vector<bool> matched(outgoing.size(),false);
for(list<pair<PDPtr,NBVertex> >::iterator it=vertices.begin();
it!=vertices.end();++it) {
// search down the tree
if(it->second.vertex) {
constructTag(channelType,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;
channelType.push_back(ix+1);
break;
}
++ix;
}
}
}
}
NBVertex::NBVertex(PrototypeVertexPtr proto) {
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;
}
}
}
}
diff --git a/Models/General/PrototypeVertex.h b/Models/General/PrototypeVertex.h
--- a/Models/General/PrototypeVertex.h
+++ b/Models/General/PrototypeVertex.h
@@ -1,474 +1,477 @@
// -*- 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 <stack>
#include "ThePEG/Helicity/Vertex/VertexBase.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Utilities/EnumIO.h"
+#include "NBodyDecayConstructorBase.fh"
//
// 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) {
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() (const pair< tPDPtr, PrototypeVertexPtr > & p1,
const pair< tPDPtr, PrototypeVertexPtr > & p2) const {
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;
/**
* Create a \f$1\to2\f$ prototype
*/
static void createPrototypes(tPDPtr inpart, VertexBasePtr vertex,
- std::stack<PrototypeVertexPtr> & prototypes);
+ std::stack<PrototypeVertexPtr> & prototypes,
+ NBodyDecayConstructorBasePtr decayCon);
/**
* Expand the prototypes by adding more legs
*/
static void expandPrototypes(PrototypeVertexPtr proto, VertexBasePtr vertex,
std::stack<PrototypeVertexPtr> & prototypes,
- const set<PDPtr> & excluded);
+ const set<PDPtr> & excluded,
+ NBodyDecayConstructorBasePtr decayCon);
/**
* 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.empty()&&y.outgoing.empty()) return true;
if(x.outgoing.size() != y.outgoing.size()) return false;
OrderedVertices::const_iterator xt = x.outgoing.begin();
OrderedVertices::const_iterator yt = y.outgoing.begin();
while(xt!=x.outgoing.end()) {
if(xt->first != yt->first) return false;
// special for identical particles
OrderedVertices::const_iterator lxt = x.outgoing.lower_bound(*xt);
OrderedVertices::const_iterator uxt = x.outgoing.upper_bound(*xt);
--uxt;
// just one particle
if(lxt==uxt) {
if(xt->second && yt->second) {
if(*(xt->second)==*(yt->second)) {
++xt;
++yt;
continue;
}
else return false;
}
else if(xt->second || yt->second)
return false;
++xt;
++yt;
}
// identical particles
else {
++uxt;
OrderedVertices::const_iterator lyt = y.outgoing.lower_bound(*xt);
OrderedVertices::const_iterator uyt = y.outgoing.upper_bound(*xt);
unsigned int nx=0;
for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {++nx;}
unsigned int ny=0;
for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {++ny;}
if(nx!=ny) return false;
vector<bool> matched(ny,false);
for(OrderedVertices::const_iterator ixt=lxt;ixt!=uxt;++ixt) {
bool found = false;
unsigned int iy=0;
for(OrderedVertices::const_iterator iyt=lyt;iyt!=uyt;++iyt) {
if(matched[iy]) {
++iy;
continue;
}
if( (!ixt->second &&!iyt->second) ||
( ixt->second&&iyt->second &&
*(ixt->second)==*(iyt->second)) ) {
matched[iy] = true;
found = true;
break;
}
++iy;
}
if(!found) return false;
}
xt=uxt;
yt=uyt;
}
}
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\to n\f$ decay
* that has been automatically generated.
*/
struct NBDiagram : public NBVertex {
/**
* Constructor taking a prototype vertex as the arguments*/
NBDiagram(PrototypeVertexPtr proto=PrototypeVertexPtr());
/**
* The type of channel
*/
vector<unsigned int> channelType;
/** 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;
};
/**
* Output operator to allow the structure to be persistently written
* @param os The output stream
* @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 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 NBDiagram
*/
inline PersistentOStream & operator<<(PersistentOStream & os,
const NBDiagram & x) {
os << x.incoming << 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 NBDiagram
*/
inline PersistentIStream & operator>>(PersistentIStream & is,
NBDiagram & x) {
is >> x.incoming >> 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 ";
for(unsigned int ix=0;ix<diag.channelType.size();++ix)
os << diag.channelType[ix] << " ";
os << "\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 */
diff --git a/Models/General/ThreeBodyDecayConstructor.cc b/Models/General/ThreeBodyDecayConstructor.cc
--- a/Models/General/ThreeBodyDecayConstructor.cc
+++ b/Models/General/ThreeBodyDecayConstructor.cc
@@ -1,363 +1,380 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ThreeBodyDecayConstructor class.
//
#include "ThreeBodyDecayConstructor.h"
#include "ThePEG/Utilities/Debug.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "Herwig++/Decay/General/GeneralThreeBodyDecayer.h"
#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
#include "Herwig++/PDT/ThreeBodyAllOnCalculator.h"
#include "ThePEG/PDT/StandardMatchers.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Utilities/Throw.h"
#include "DecayConstructor.h"
#include "WeakCurrentDecayConstructor.h"
using namespace Herwig;
IBPtr ThreeBodyDecayConstructor::clone() const {
return new_ptr(*this);
}
IBPtr ThreeBodyDecayConstructor::fullclone() const {
return new_ptr(*this);
}
void ThreeBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
os << interOpt_ << widthOpt_ << intOpt_ << relErr_
<< includeIntermediatePhotons_;
}
void ThreeBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
is >> interOpt_ >> widthOpt_ >> intOpt_ >> relErr_
>> includeIntermediatePhotons_;
}
ClassDescription<ThreeBodyDecayConstructor>
ThreeBodyDecayConstructor::initThreeBodyDecayConstructor;
// Definition of the static class description member.
void ThreeBodyDecayConstructor::Init() {
static ClassDocumentation<ThreeBodyDecayConstructor> documentation
("The ThreeBodyDecayConstructor class constructs the three body decay modes");
static Switch<ThreeBodyDecayConstructor,bool> interfaceIncludeIntermediatePhotons
("IncludeIntermediatePhotons",
"Whether or not ot allow intermediate photons",
&ThreeBodyDecayConstructor::includeIntermediatePhotons_, false, false, false);
static SwitchOption interfaceIncludeIntermediatePhotonsYes
(interfaceIncludeIntermediatePhotons,
"Yes",
"Include them",
true);
static SwitchOption interfaceIncludeIntermediatePhotonsNo
(interfaceIncludeIntermediatePhotons,
"No",
"Don't include them",
false);
static Switch<ThreeBodyDecayConstructor,unsigned int> interfaceWidthOption
("WidthOption",
"Option for the treatment of the widths of the intermediates",
&ThreeBodyDecayConstructor::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<ThreeBodyDecayConstructor,unsigned int> interfaceIntermediateOption
("IntermediateOption",
"Option for the inclusion of intermediates in the event",
&ThreeBodyDecayConstructor::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);
static Switch<ThreeBodyDecayConstructor,unsigned int> interfaceIntegrationOption
("IntegrationOption",
"Option for the integration of the partial width",
&ThreeBodyDecayConstructor::intOpt_, 1, false, false);
static SwitchOption interfaceIntegrationOptionAllPoles
(interfaceIntegrationOption,
"AllPoles",
"Include all potential poles",
0);
static SwitchOption interfaceIntegrationOptionShallowestPole
(interfaceIntegrationOption,
"ShallowestPole",
"Only include the shallowest pole",
1);
static Parameter<ThreeBodyDecayConstructor,double> interfaceRelativeError
("RelativeError",
"The relative error for the GQ integration",
&ThreeBodyDecayConstructor::relErr_, 1e-2, 1e-10, 1.,
false, false, Interface::limited);
}
void ThreeBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
if( particles.empty() ) return;
// special for weak decays
for(unsigned int ix=0;ix<decayConstructor()->decayConstructors().size();++ix) {
Ptr<Herwig::WeakCurrentDecayConstructor>::pointer
weak = dynamic_ptr_cast<Ptr<Herwig::WeakCurrentDecayConstructor>::pointer >
(decayConstructor()->decayConstructors()[ix]);
if(!weak) continue;
weakMassCut_ = max(weakMassCut_,weak->massCut());
}
NBodyDecayConstructorBase::DecayList(particles);
}
GeneralThreeBodyDecayerPtr ThreeBodyDecayConstructor::
createDecayer(vector<TBDiagram> & diagrams, bool inter,
double symfac) const {
if(diagrams.empty()) return GeneralThreeBodyDecayerPtr();
// extract the external particles for the process
PDPtr incoming = getParticleData(diagrams[0].incoming);
// outgoing particles
OrderedParticles outgoing;
outgoing.insert(getParticleData(diagrams[0].outgoing ));
outgoing.insert(getParticleData(diagrams[0].outgoingPair.first ));
outgoing.insert(getParticleData(diagrams[0].outgoingPair.second));
// get the object name
string objectname ("/Herwig/Decays/");
string classname = DecayerClassName(incoming, outgoing, objectname);
if(classname=="") return GeneralThreeBodyDecayerPtr();
// create the object
GeneralThreeBodyDecayerPtr decayer =
dynamic_ptr_cast<GeneralThreeBodyDecayerPtr>
(generator()->preinitCreate(classname, objectname));
// set up the decayer and return if doesn't work
vector<PDPtr> outVector(outgoing.begin(),outgoing.end());
if(!decayer->setDecayInfo(incoming,outVector,diagrams,symfac))
return GeneralThreeBodyDecayerPtr();
// set decayer options from base class
setDecayerInterfaces(objectname);
// options for partial width integration
ostringstream value;
value << intOpt_;
generator()->preinitInterface(objectname, "PartialWidthIntegration", "set",
value.str());
value.str("");
value << relErr_;
generator()->preinitInterface(objectname, "RelativeError", "set",
value.str());
// set the width option
value.str("");
value << widthOpt_;
generator()->preinitInterface(objectname, "WidthOption", "set", value.str());
// set the intermediates option
value.str("");
value << inter;
generator()->preinitInterface(objectname, "GenerateIntermediates", "set",
value.str());
// initialize the decayer
decayer->init();
// return the decayer
return decayer;
}
string ThreeBodyDecayConstructor::
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(ns==1&&nf==2) classname += "StoSFFDecayer";
else if(nf==2&&nv==1) classname += "StoFFVDecayer";
else classname = "";
}
else if(incoming->iSpin()==PDT::Spin1Half) {
if(nf==3) classname += "FtoFFFDecayer";
else if(nf==1&&nv==2) classname += "FtoFVVDecayer";
else classname = "";
}
else if(incoming->iSpin()==PDT::Spin1) {
if(nf==2&&nv==1) classname += "VtoFFVDecayer";
else classname = "";
}
else {
classname="";
}
return classname;
}
void ThreeBodyDecayConstructor::
createDecayMode(vector<NBDiagram> & mode,
bool possibleOnShell,
double symfac) {
// convert the diagrams from the general to the three body structure
vector<TBDiagram> diagrams;
for(unsigned int iy=0;iy<mode.size();++iy) {
diagrams.push_back(TBDiagram(mode[iy]));
// determine the type
diagrams.back().channelType = TBDiagram::Channel(mode[iy].channelType[0]-1);
// remove weak processes simulated using the weak current
if(weakMassCut_>ZERO && diagrams.back().intermediate &&
abs(diagrams.back().intermediate->id())==ParticleID::Wplus) {
Energy deltaM =
getParticleData(diagrams.back().incoming)->mass() -
getParticleData(diagrams.back().outgoing)->mass();
if(deltaM<weakMassCut_) diagrams.pop_back();
}
// remove intermediate photons
else if(!includeIntermediatePhotons_ &&
diagrams.back().intermediate &&
abs(diagrams.back().intermediate->id())==ParticleID::gamma)
diagrams.pop_back();
}
if(diagrams.empty()) return;
// check if possible on-shell internal particles
bool inter = interOpt_ == 1 || (interOpt_ == 0 && possibleOnShell);
// incoming particle
tPDPtr inpart = getParticleData(diagrams[0].incoming);
// outgoing particles
OrderedParticles outgoing;
outgoing.insert(getParticleData(diagrams[0].outgoing));
outgoing.insert(getParticleData(diagrams[0].outgoingPair.first ));
outgoing.insert(getParticleData(diagrams[0].outgoingPair.second));
// sort out ordering and labeling of diagrams
vector<PDPtr> outVector(outgoing.begin(),outgoing.end());
for(unsigned int ix=0;ix<diagrams.size();++ix) {
if( ( diagrams[ix].channelType == TBDiagram::channel23 &&
outVector[1]->id() != diagrams[ix].outgoingPair.first) ||
( diagrams[ix].channelType == TBDiagram::channel13 &&
outVector[0]->id() != diagrams[ix].outgoingPair.first) ||
( diagrams[ix].channelType == TBDiagram::channel12 &&
outVector[0]->id() != diagrams[ix].outgoingPair.first) )
swap(diagrams[ix].outgoingPair.first, diagrams[ix].outgoingPair.second);
}
// construct the tag for the decay mode
string tag = inpart->name() + "->";
unsigned int iprod=0;
for(OrderedParticles::const_iterator it = outgoing.begin();
it != outgoing.end(); ++it) {
++iprod;
tag += (**it).name();
if(iprod != 3) tag += ",";
}
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;
}
// create mode if needed
if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
// create the decayer
GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
if(!decayer) {
if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
<< tag << " so mode not created\n";
return;
}
tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
if(ndm) {
generator()->preinitInterface(ndm, "Decayer", "set",
decayer->fullName());
generator()->preinitInterface(ndm, "OnOff", "set", "On");
if(!ndm->decayer()) {
generator()->log() << "Can't set the decayer for "
<< tag << " so mode not created \n";
return;
}
OrderedParticles::const_iterator pit=outgoing.begin();
tPDPtr pa = *pit; ++pit;
tPDPtr pb = *pit; ++pit;
tPDPtr pc = *pit;
Energy width =
decayer->partialWidth(make_pair(inpart,inpart->mass()),
make_pair(pa,pa->mass()) ,
make_pair(pb,pb->mass()) ,
make_pair(pc,pc->mass()));
if ( Debug::level > 1 )
generator()->log() << "Partial width is: " << width / GeV << " GeV\n";
setBranchingRatio(ndm, width);
if(ndm->brat()<decayConstructor()->minimumBR()) {
generator()->preinitInterface(decayer->fullName(),
"Initialize", "set","0");
}
// incoming particle is now unstable
inpart->stable(false);
if(Debug::level > 1 ) {
generator()->log() << "Calculated partial width for mode "
<< tag << " is " << width/GeV << "\n";
}
}
else
throw NBodyDecayConstructorError()
<< "ThreeBodyDecayConstructor::createDecayMode - Needed to create "
<< "new decaymode but one could not be created for the tag "
<< tag << Exception::warning;
}
else if( dm ) {
if(dm->brat()<decayConstructor()->minimumBR()) {
return;
}
if((dm->decayer()->fullName()).find("Mambo") != string::npos) {
// create the decayer
GeneralThreeBodyDecayerPtr decayer = createDecayer(diagrams,inter,symfac);
- if(!decayer) {
- if(Debug::level > 1 ) generator()->log() << "Can't create the decayer for "
- << tag << " so mode not created\n";
- return;
+ if(decayer) {
+ generator()->preinitInterface(dm, "Decayer", "set",
+ decayer->fullName());
+ // check not zero
+ OrderedParticles::const_iterator pit=outgoing.begin();
+ tPDPtr pa = *pit; ++pit;
+ tPDPtr pb = *pit; ++pit;
+ tPDPtr pc = *pit;
+ Energy width =
+ decayer->partialWidth(make_pair(inpart,inpart->mass()),
+ make_pair(pa,pa->mass()) ,
+ make_pair(pb,pb->mass()) ,
+ make_pair(pc,pc->mass()));
+ if(width/(dm->brat()*inpart->width())<1e-10) {
+ string message = "Herwig++ calculation of the partial width for the decay mode "
+ + inpart->PDGName() + " -> " + pa->PDGName() + " " + pb->PDGName() + " " + pc->PDGName()
+ + " is zero.\n This will cause problems with the calculation of"
+ + " spin correlations.\n It is probably due to inconsistent parameters"
+ + " and decay modes being passed to Herwig++ via the SLHA file.\n"
+ + " Zeroing the branching ratio for this mode.";
+ setBranchingRatio(dm,ZERO);
+ generator()->logWarning(NBodyDecayConstructorError(message,Exception::warning));
+ }
}
- generator()->preinitInterface(dm, "Decayer", "set",
- decayer->fullName());
// incoming particle is now unstable
inpart->stable(false);
}
}
//update CC mode if it exists
if( inpart->CC() ) inpart->CC()->synchronize();
}
diff --git a/Models/General/TwoBodyDecayConstructor.cc b/Models/General/TwoBodyDecayConstructor.cc
--- a/Models/General/TwoBodyDecayConstructor.cc
+++ b/Models/General/TwoBodyDecayConstructor.cc
@@ -1,365 +1,379 @@
// -*- C++ -*-
//
// TwoBodyDecayConstructor.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TwoBodyDecayConstructor class.
//
#include "TwoBodyDecayConstructor.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include "Herwig++/Decay/General/GeneralTwoBodyDecayer.h"
#include "Herwig++/Models/StandardModel/StandardModel.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "DecayConstructor.h"
#include "ThePEG/Utilities/Throw.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractFFTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSTVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractSSSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFSVertex.fh"
#include "ThePEG/Helicity/Vertex/AbstractRFVVertex.fh"
using namespace Herwig;
using ThePEG::Helicity::VertexBasePtr;
IBPtr TwoBodyDecayConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoBodyDecayConstructor::fullclone() const {
return new_ptr(*this);
}
NoPIOClassDescription<TwoBodyDecayConstructor>
TwoBodyDecayConstructor::initTwoBodyDecayConstructor;
// Definition of the static class description member.
void TwoBodyDecayConstructor::Init() {
static ClassDocumentation<TwoBodyDecayConstructor> documentation
("The TwoBodyDecayConstructor implements to creation of 2 body decaymodes "
"and decayers that do not already exist for the given set of vertices.");
}
void TwoBodyDecayConstructor::DecayList(const set<PDPtr> & particles) {
if( particles.empty() ) return;
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
unsigned int nv(model->numberOfVertices());
for(set<PDPtr>::const_iterator ip=particles.begin();
ip!=particles.end();++ip) {
tPDPtr parent = *ip;
if ( Debug::level > 0 )
Repository::cout() << "Constructing 2-body decays for "
<< parent->PDGName() << '\n';
for(unsigned int iv = 0; iv < nv; ++iv) {
if(excluded(model->vertex(iv)) ||
model->vertex(iv)->getNpoint()>3) continue;
for(unsigned int il = 0; il < 3; ++il) {
set<TwoBodyDecay> decays =
createModes(parent, model->vertex(iv), il);
if( !decays.empty() ) createDecayMode(decays);
}
}
}
}
set<TwoBodyDecay> TwoBodyDecayConstructor::
createModes(tPDPtr inpart, VertexBasePtr vertex,
unsigned int list) {
if( !vertex->isIncoming(inpart) || vertex->getNpoint() != 3 )
return set<TwoBodyDecay>();
Energy m1(inpart->mass());
tPDPtr ccpart = inpart->CC() ? inpart->CC() : inpart;
long id = ccpart->id();
tPDVector decaylist = vertex->search(list, ccpart);
set<TwoBodyDecay> decays;
tPDVector::size_type nd = decaylist.size();
for( tPDVector::size_type i = 0; i < nd; i += 3 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]);
if( pb->id() == id ) swap(pa, pb);
if( pc->id() == id ) swap(pa, pc);
//allowed on-shell decay?
if( m1 <= pb->mass() + pc->mass() ) continue;
//vertices are defined with all particles incoming
decays.insert( TwoBodyDecay(inpart,pb, pc, vertex) );
}
return decays;
}
GeneralTwoBodyDecayerPtr TwoBodyDecayConstructor::createDecayer(TwoBodyDecay decay) {
string name;
using namespace Helicity::VertexType;
PDT::Spin in = decay.parent_->iSpin();
// PDT::Spin out1 = decay.children_.first ->iSpin();
PDT::Spin out2 = decay.children_.second->iSpin();
switch(decay.vertex_->getName()) {
case FFV :
if(in == PDT::Spin1Half) {
name = "FFVDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "VFFDecayer";
}
break;
case FFS :
if(in == PDT::Spin1Half) {
name = "FFSDecayer";
if(out2==PDT::Spin1Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SFFDecayer";
}
break;
case VVS :
if(in == PDT::Spin1) {
name = "VVSDecayer";
if(out2==PDT::Spin1)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "SVVDecayer";
}
break;
case VSS :
if(in == PDT::Spin1) {
name = "VSSDecayer";
}
else {
name = "SSVDecayer";
if(out2==PDT::Spin0)
swap(decay.children_.first,decay.children_.second);
}
break;
case VVT :
name = in==PDT::Spin2 ? "TVVDecayer" : "Unknown";
break;
case FFT :
name = in==PDT::Spin2 ? "TFFDecayer" : "Unknown";
break;
case SST :
name = in==PDT::Spin2 ? "TSSDecayer" : "Unknown";
break;
case SSS :
name = "SSSDecayer";
break;
case VVV :
name = "VVVDecayer";
break;
case RFS :
if(in==PDT::Spin1Half) {
name = "FRSDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else if(in==PDT::Spin0) {
name = "SRFDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else {
name = "Unknown";
}
break;
case RFV :
if(in==PDT::Spin1Half) {
name = "FRVDecayer";
if(out2==PDT::Spin3Half)
swap(decay.children_.first,decay.children_.second);
}
else
name = "Unknown";
break;
default : Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
}
if(name=="Unknown")
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
ostringstream fullname;
fullname << "/Herwig/Decays/" << name << "_" << decay.parent_->PDGName()
<< "_" << decay.children_.first ->PDGName()
<< "_" << decay.children_.second->PDGName();
string classname = "Herwig::" + name;
GeneralTwoBodyDecayerPtr decayer;
decayer = dynamic_ptr_cast<GeneralTwoBodyDecayerPtr>
(generator()->preinitCreate(classname,fullname.str()));
if(!decayer)
Throw<NBodyDecayConstructorError>()
<< "Error: Cannot assign " << decay.vertex_->fullName() << " to a decayer. "
<< "Decay is " << decay.parent_->PDGName() << " -> "
<< decay.children_.first ->PDGName() << " "
<< decay.children_.second->PDGName() << Exception::runerror;
// set the strong coupling for radiation
generator()->preinitInterface(decayer, "Coupling", "set", showerAlpha_);
// get the vertices for radiation from the external legs
VertexBasePtr inRad = radiationVertex(decay.parent_);
vector<VertexBasePtr> outRad;
outRad.push_back(radiationVertex(decay.children_.first ));
outRad.push_back(radiationVertex(decay.children_.second));
// get any contributing 4 point vertices
VertexBasePtr fourRad = radiationVertex(decay.parent_, decay.children_);
// set info on decay
decayer->setDecayInfo(decay.parent_,decay.children_,decay.vertex_,
inRad,outRad,fourRad);
// initialised the decayer
setDecayerInterfaces(fullname.str());
decayer->init();
return decayer;
}
void TwoBodyDecayConstructor::
createDecayMode(set<TwoBodyDecay> & decays) {
tPDPtr inpart = decays.begin()->parent_;
set<TwoBodyDecay>::iterator dend = decays.end();
for( set<TwoBodyDecay>::iterator dit = decays.begin();
dit != dend; ++dit ) {
tPDPtr pb((*dit).children_.first), pc((*dit).children_.second);
string tag = inpart->name() + "->" + pb->name() + "," +
pc->name() + ";";
// Does it exist already ?
tDMPtr dm = generator()->findDecayMode(tag);
// Check if tag is one that should be disabled
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");
}
continue;
}
// now create DecayMode objects that do not already exist
if( createDecayModes() && (!dm || inpart->id() == ParticleID::h0) ) {
tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
if(ndm) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(*dit);
if(!decayer) continue;
generator()->preinitInterface(ndm, "Decayer", "set",
decayer->fullName());
generator()->preinitInterface(ndm, "OnOff", "set", "On");
Energy width =
decayer->partialWidth(make_pair(inpart,inpart->mass()),
make_pair(pb,pb->mass()) ,
make_pair(pc,pc->mass()));
setBranchingRatio(ndm, width);
if(ndm->brat()<decayConstructor()->minimumBR()) {
generator()->preinitInterface(decayer->fullName(),
"Initialize", "set","0");
}
}
else
Throw<NBodyDecayConstructorError>()
<< "TwoBodyDecayConstructor::createDecayMode - Needed to create "
<< "new decaymode but one could not be created for the tag "
<< tag << Exception::warning;
}
else if( dm ) {
if(dm->brat()<decayConstructor()->minimumBR()) {
continue;
}
if((dm->decayer()->fullName()).find("Mambo") != string::npos) {
inpart->stable(false);
GeneralTwoBodyDecayerPtr decayer=createDecayer(*dit);
if(!decayer) continue;
generator()->preinitInterface(dm, "Decayer", "set",
decayer->fullName());
+ Energy width =
+ decayer->partialWidth(make_pair(inpart,inpart->mass()),
+ make_pair(pb,pb->mass()) ,
+ make_pair(pc,pc->mass()));
+ if(width/(dm->brat()*inpart->width())<1e-10) {
+ string message = "Herwig++ calculation of the partial width for the decay mode "
+ + inpart->PDGName() + " -> " + pb->PDGName() + " " + pc->PDGName()
+ + " is zero.\n This will cause problems with the calculation of"
+ + " spin correlations.\n It is probably due to inconsistent parameters"
+ + " and decay modes being passed to Herwig++ via the SLHA file.\n"
+ + " Zeroing the branching ratio for this mode.";
+ setBranchingRatio(dm,ZERO);
+ generator()->logWarning(NBodyDecayConstructorError(message,Exception::warning));
+ }
}
}
}
// update CC mode if it exists
if( inpart->CC() ) inpart->CC()->synchronize();
}
VertexBasePtr TwoBodyDecayConstructor::radiationVertex(tPDPtr particle, tPDPair children) {
tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
map<tPDPtr,VertexBasePtr>::iterator rit = radiationVertices_.find(particle);
tPDPtr cc = particle->CC() ? particle->CC() : particle;
if(children==tPDPair() && rit!=radiationVertices_.end()) return rit->second;
unsigned int nv(model->numberOfVertices());
tPDPtr gluon = getParticleData(ParticleID::g);
// look for radiation vertices for incoming and outgoing particles
for(unsigned int iv=0;iv<nv;++iv) {
VertexBasePtr vertex = model->vertex(iv);
// look for 3 point vertices
if (children==tPDPair()){
if( !vertex->isIncoming(particle) || vertex->getNpoint() != 3 ||
!vertex->isOutgoing(particle) || !vertex->isOutgoing(gluon)) continue;
for(unsigned int list=0;list<3;++list) {
tPDVector decaylist = vertex->search(list, particle);
for( tPDVector::size_type i = 0; i < decaylist.size(); i += 3 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i + 1]), pc(decaylist[i + 2]);
if( pb->id() == ParticleID::g ) swap(pa, pb);
if( pc->id() == ParticleID::g ) swap(pa, pc);
if( pb->id() != particle->id()) swap(pb, pc);
if( pa->id() != ParticleID::g) continue;
if( pb != particle) continue;
if( pc != cc) continue;
radiationVertices_[particle] = vertex;
return vertex;
}
}
}
// look for 4 point vertex including a gluon
else {
if( !vertex->isIncoming(particle) || vertex->getNpoint()!=4 ||
!vertex->isOutgoing(children.first) || !vertex->isOutgoing(children.second) ||
!vertex->isOutgoing(gluon)) continue;
for(unsigned int list=0;list<4;++list) {
tPDVector decaylist = vertex->search(list, particle);
for( tPDVector::size_type i = 0; i < decaylist.size(); i += 4 ) {
tPDPtr pa(decaylist[i]), pb(decaylist[i+1]), pc(decaylist[i+2]), pd(decaylist[i+3]);
// order so that a = g, b = parent
if( pb->id() == ParticleID::g ) swap(pa, pb);
if( pc->id() == ParticleID::g ) swap(pa, pc);
if( pd->id() == ParticleID::g ) swap(pa, pd);
if( pc->id() == particle->id()) swap(pb, pc);
if( pd->id() == particle->id()) swap(pb, pd);
if( pa->id() != ParticleID::g) continue;
if( pb->id() != particle->id()) continue;
if( !((abs(pd->id()) == abs(children. first->id()) &&
abs(pc->id()) == abs(children.second->id())) ||
(abs(pc->id()) == abs(children. first->id()) &&
abs(pd->id()) == abs(children.second->id()))))
continue;
return vertex;
}
}
}
}
return VertexBasePtr();
}
diff --git a/Models/LH/LHModel.h b/Models/LH/LHModel.h
--- a/Models/LH/LHModel.h
+++ b/Models/LH/LHModel.h
@@ -1,289 +1,289 @@
// -*- C++ -*-
#ifndef HERWIG_LHModel_H
#define HERWIG_LHModel_H
//
// This is the declaration of the LHModel class.
//
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSSVertex.h"
#include "LHModel.fh"
namespace Herwig {
using namespace ThePEG;
/**
* The LHModel class is the main class for the implementation of the Little Higgs model in Herwig++.
* In inherits from the StandardModel class and implements the calcuation of the couplings
* and masses in the Little Higgs model and storage of the additional couplings.
*
* @see \ref LHModelInterfaces "The interfaces"
* defined for LHModel.
*/
-class LHModel: public StandardModel {
+class LHModel: public BSMModel {
public:
/**
* The default constructor.
*/
LHModel();
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:
/**
* Access to the parameters of the model
*/
//@{
/**
* The \f$\lambda_1\f$ top Yukawa coupling
*/
double lambda1() const {return _lambda1;}
/**
* The \f$\lambda_2\f$ top Yukawa coupling
*/
double lambda2() const {return _lambda2;}
/**
* The sine of the \f$\theta\f$ mixing angle
*/
double sinTheta() const {return _s;}
/**
* The cosine of the \f$\theta\f$ mixing angle
*/
double cosTheta() const {return _c;}
/**
* The sine of the \f$\theta'\f$ mixing angle
*/
double sinThetaPrime() const {return _sp;}
/**
* The cosine of the \f$\theta'\f$ mixing angle
*/
double cosThetaPrime() const {return _cp;}
/**
* The sine of the Higgs mixing angle
*/
double sinTheta0() const {return _s0;}
/**
* The cosine of the Higgs mixing angle
*/
double cosTheta0() const {return _c0;}
/**
* The sine of the pseudoscalar Higgs mixing angle
*/
double sinThetaP() const {return _sP;}
/**
* The cosine of the pseudoscalar Higgs mixing angle
*/
double cosThetaP() const {return _cP;}
/**
* The sine of the charged Higgs mixing angle
*/
double sinThetaPlus() const {return _sPlus;}
/**
* The cosine of the charged Higgs mixing angle
*/
double cosThetaPlus() const {return _cPlus;}
/**
* The vacuum expection value
*/
Energy vev() const {return _v;}
/**
* The vacuum expection value
*/
Energy vevPrime() const {return _v*_vacratio;}
/**
* The \f$f\f$ scale of the non-linear \f$\sigma\f$-model
*/
Energy f() const {return _f;}
//@}
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
LHModel & operator=(const LHModel &);
private:
/**
* Parameters for the model
*/
//@{
/**
* The value of \f$\cot\theta\f$ for the mixing with the \f$g\f$ coupling
*/
double _cott;
/**
* The value of \f$\tan\theta'\f$ for the mixing with the \f$g'\f$ coupling
*/
double _tantp;
/**
* The vacuum expection valve
*/
Energy _v;
/**
* The ratio \f$\lambda_2/\lambda_1\f$
*/
double _lamratio;
/**
* The mass of the lightest Higgs boson
*/
Energy _mH;
/**
* The ratio of the vacuum exception values \f$v'/v\f$
*/
double _vacratio;
/**
* The scale for the non-linear \f$\sigma\f$-model
*/
Energy _f;
/**
* The top Yukawa couplings
*/
double _lambda1,_lambda2;
/**
* The sine of the \f$\theta\f$ mixing angle
*/
double _s;
/**
* The cosine of the \f$\theta\f$ mixing angle
*/
double _c;
/**
* The sine of the \f$\theta'\f$ mixing angle
*/
double _sp;
/**
* The cosine of the \f$\theta'\f$ mixing angle
*/
double _cp;
/**
* The sine of the Higgs mixing angle
*/
double _s0;
/**
* The cosine of the Higgs mixing angle
*/
double _c0;
/**
* The sine of the pseudoscalar Higgs mixing angle
*/
double _sP;
/**
* The cosine of the pseudoscalar Higgs mixing angle
*/
double _cP;
/**
* The sine of the charged Higgs mixing angle
*/
double _sPlus;
/**
* The cosine of the charged Higgs mixing angle
*/
double _cPlus;
//@}
/**
* Additional vertices
*/
//@{
/**
* WHH Vertex
*/
AbstractVSSVertexPtr WHHVertex_;
//@}
};
}
#endif /* HERWIG_LHModel_H */
diff --git a/Models/LHTP/LHTPModel.h b/Models/LHTP/LHTPModel.h
--- a/Models/LHTP/LHTPModel.h
+++ b/Models/LHTP/LHTPModel.h
@@ -1,302 +1,302 @@
// -*- C++ -*-
#ifndef THEPEG_LHTPModel_H
#define THEPEG_LHTPModel_H
//
// This is the declaration of the LHTPModel class.
//
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "LHTPModel.fh"
namespace Herwig {
/**
* The LHTPModel class is the main class for the
* implementation of the Little Higgs model with T-parity
*
* @see \ref LHTPModelInterfaces "The interfaces"
* defined for LHTPModel.
*/
-class LHTPModel: public StandardModel {
+class LHTPModel: public BSMModel {
public:
/**
* The default constructor.
*/
LHTPModel();
/**
* Access to the parameters of the model
*/
//@{
/**
* The vacuum expection value
*/
Energy vev() const { return v_; }
/**
* The \f$f\f$ scale of the non-linear \f$\sigma\f$-model
*/
Energy f() const { return f_; }
/**
* \f$\sin\alpha\f$
*/
double sinAlpha() const { return salpha_; }
/**
* \f$\cos\alpha\f$
*/
double cosAlpha() const { return calpha_; }
/**
* \f$\sin\beta\f$
*/
double sinBeta() const { return sbeta_; }
/**
* \f$\cos\beta\f$
*/
double cosBeta() const { return cbeta_; }
/**
* \f$\sin\theta_H\f$
*/
double sinThetaH() const { return sthetaH_; }
/**
* \f$\cos\theta_H\f$
*/
double cosThetaH() const { return cthetaH_; }
/**
* \f$\sin\theta_L\f$
*/
double sinThetaL() const { return sL_;}
/**
* \f$\cos\theta_L\f$
*/
double cosThetaL() const { return cL_;}
/**
* \f$\sin\theta_R\f$
*/
double sinThetaR() const { return sR_;}
/**
* \f$\cos\theta_R\f$
*/
double cosThetaR() const { return cR_;}
//@}
/**
* Yukwawa for T-odd fermions
*/
//@{
/**
* The \f$\kappa_q\f$ parameter which controls the properties of the
* T-odd quarks
*/
double kappaQuark() const { return kappaQuark_;}
/**
* The \f$\kappa_\ell\f$ parameter which controls the properties of the
* T-odd leptons
*/
double kappaLepton() const { return kappaLepton_;}
//@}
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:
/**
* Calculate the mixing in the top sector of the model
* and the masses of the T-odd and T-even heavy tops
* The mixings are calculated by solving Eqns 2.22 and 2.24 of hep-ph/0506042
* for \f$\lambda_1$ and \f$\lambda_2\f$ given the input value of \f$\sin\alpha\f$
* and the top mass.
*/
void topMixing(Energy & MTp, Energy & MTm);
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.
*/
LHTPModel & operator=(const LHTPModel &);
private:
/**
* The constant for the non-linear \f$\sigma\f$ model
*/
Energy f_;
/**
* @name The mixing in the top quark sector
*/
//@{
/**
* \f$\sin\alpha\f$, taken as an input
*/
double salpha_;
/**
* \f$\cos\alpha\f$
*/
double calpha_;
/**
* \f$\sin\beta\f$
*/
double sbeta_;
/**
* \f$\cos\beta\f$
*/
double cbeta_;
//@}
/**
* @name Mixing of the heavy photon and Z
*/
//@{
/**
* \f$\sin\theta_H\f$
*/
double sthetaH_;
/**
* \f$\cos\theta_H\f$
*/
double cthetaH_;
//@}
/**
* @name Mixings in the top sector
*/
//@{
/**
* \f$\sin\theta_L\f$
*/
double sL_;
/**
* \f$\cos\theta_L\f$
*/
double cL_;
/**
* \f$\sin\theta_R\f$
*/
double sR_;
/**
* \f$\cos\theta_R\f$
*/
double cR_;
//@}
/**
* The \f$\kappa_q\f$ parameter which controls the properties of the
* T-odd quarks
*/
double kappaQuark_;
/**
* The \f$\kappa_\ell\f$ parameter which controls the properties of the
* T-odd leptons
*/
double kappaLepton_;
/**
* The mass of the Standard Model higgs
*/
Energy mh_;
/**
* The vacuum expection valve
*/
Energy v_;
/**
* The \f$g\f$ coupling
*/
double g_;
/**
* the \f$g'\f$ coupling
*/
double gp_;
/**
* Method for evaluating the masses
*/
bool approximate_;
/**
* WHH Vertex
*/
AbstractVSSVertexPtr WHHVertex_;
};
}
#endif /* THEPEG_LHTPModel_H */
diff --git a/Models/Sextet/SextetModel.h b/Models/Sextet/SextetModel.h
--- a/Models/Sextet/SextetModel.h
+++ b/Models/Sextet/SextetModel.h
@@ -1,308 +1,308 @@
// -*- C++ -*-
#ifndef HERWIG_SextetModel_H
#define HERWIG_SextetModel_H
//
// This is the declaration of the SextetModel class.
//
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVVVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractFFSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVSSVertex.h"
#include "ThePEG/Helicity/Vertex/AbstractVVSSVertex.h"
#include "SextetModel.fh"
namespace Herwig {
using namespace ThePEG;
/** \ingroup Models
*
* This class is used instead of the StandardModel class for the
*
*
* @see \ref SextetModelInterfaces "The interfaces"
* defined for SextetModel.
*/
-class SextetModel: public StandardModel {
+class SextetModel: public BSMModel {
public:
/**
* The default constructor.
*/
SextetModel() : g1L_(3,0.), g1R_(3,0.), g1pR_(3,0.), g1ppR_(3,0.),
g2_(3,0.), g2p_(3,0.), g3L_(3,0.),
enableScalarSingletY43_(false),enableScalarSingletY13_(false),
enableScalarSingletY23_(false),enableScalarTripletY13_(false),
enableVectorDoubletY16_(false),enableVectorDoubletY56_(false) {
useMe();
}
/**
* Access to the couplings
*/
//@{
/**
* The \f$SU(2)\f$ quark-doublet coupling to \f$\Phi_{6,1,1/3}\f$
*/
const vector<double> & g1L() const {return g1L_;}
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,1/3}\f$
*/
const vector<double> & g1R() const {return g1R_;}
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,-2/3}\f$
*/
const vector<double> & g1pR() const {return g1pR_;}
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,4/3}\f$
*/
const vector<double> & g1ppR() const {return g1ppR_;}
/**
* The coupling to \f$V^\mu_{6,2,-1/6}\f$
*/
const vector<double> & g2() const {return g2_;}
/**
* The coupling to \f$V^\mu_{6,2,5/6}\f$
*/
const vector<double> & g2p() const {return g2p_;}
/**
* Coupling to \f$\Phi_{6,3,1/3}\f$
*/
const vector<double> & g3L() const {return g3L_;}
//@}
/**
* Switches to decide which particles to include
*/
//@{
/**
* Scalar Singlet \f$Y = 4/3\f$
*/
bool ScalarSingletY43Enabled() const {return enableScalarSingletY43_;}
/**
* Scalar Singlet \f$Y = -1/3\f$
*/
bool ScalarSingletY13Enabled() const {return enableScalarSingletY13_;}
/**
* Scalar Singlet \f$Y = -2/3\f$
*/
bool ScalarSingletY23Enabled() const {return enableScalarSingletY23_;}
/**
* Scalar Triplet \f$Y = 1/3\f$
*/
bool ScalarTripletY13Enabled() const {return enableScalarTripletY13_;}
/**
* Vector Doublet \f$Y = -1/6\f$
*/
bool VectorDoubletY16Enabled() const {return enableVectorDoubletY16_;}
/**
* Vector Doublet \f$Y = 5/6\f$
*/
bool VectorDoubletY56Enabled() const {return enableVectorDoubletY56_;}
//@}
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();
/**
* Member to implement the command to enable particular diquarks
*/
string doEnable(string command);
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SextetModel & operator=(const SextetModel &);
private:
/**
* Pointers to the vertex objects
*/
//@{
/**
* Pointer to the object handling the strong coupling of a
* vector sextet to one gluon
*/
AbstractVVVVertexPtr VVVVertex_;
/**
* Pointer to the object handling the strong coupling of a
* vector sextet to two gluons
*/
AbstractVVVVVertexPtr VVVVVertex_;
/**
* Pointer to the object handling the strong coupling of a
* scalar sextet to one gluon
*/
AbstractVSSVertexPtr VSSVertex_;
/**
* Pointer to the object handling the strong coupling of a
* scalar sextet to two gluons
*/
AbstractVVSSVertexPtr VVSSVertex_;
/**
* Pointer to the object handling the coupling of two quarks
* to a vector sextet
*/
AbstractFFVVertexPtr FFVVertex_;
/**
* Pointer to the object handling the coupling of two quarks
* to a scalar sextet
*/
AbstractFFSVertexPtr FFSVertex_;
//@}
/**
* Couplings
*/
//@{
/**
* The \f$SU(2)\f$ quark-doublet coupling to \f$\Phi_{6,1,1/3}\f$
*/
vector<double> g1L_;
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,1/3}\f$
*/
vector<double> g1R_;
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,-2/3}\f$
*/
vector<double> g1pR_;
/**
* The \f$SU(2)\f$ singlet coupling to \f$\Phi_{6,1,4/3}\f$
*/
vector<double> g1ppR_;
/**
* The coupling to \f$V^\mu_{6,2,-1/6}\f$
*/
vector<double> g2_;
/**
* The coupling to \f$V^\mu_{6,2,5/6}\f$
*/
vector<double> g2p_;
/**
* Coupling to \f$\Phi_{6,3,1/3}\f$
*/
vector<double> g3L_;
//@}
/**
* Switches to decide which particles to include
*/
//@{
/**
* Scalar Singlet \f$Y = 4/3\f$
*/
bool enableScalarSingletY43_;
/**
* Scalar Singlet \f$Y = -1/3\f$
*/
bool enableScalarSingletY13_;
/**
* Scalar Singlet \f$Y = -2/3\f$
*/
bool enableScalarSingletY23_;
/**
* Scalar Triplet \f$Y = 1/3\f$
*/
bool enableScalarTripletY13_;
/**
* Vector Doublet \f$Y = -1/6\f$
*/
bool enableVectorDoubletY16_;
/**
* Vector Doublet \f$Y = 5/6\f$
*/
bool enableVectorDoubletY56_;
//@}
};
}
#endif /* HERWIG_SextetModel_H */
diff --git a/Models/Susy/SSGNGVertex.cc b/Models/Susy/SSGNGVertex.cc
--- a/Models/Susy/SSGNGVertex.cc
+++ b/Models/Susy/SSGNGVertex.cc
@@ -1,293 +1,304 @@
// -*- C++ -*-
//
// SSGNGVertex.cc is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the SSGNGVertex class.
//
#include "SSGNGVertex.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "Herwig++/Models/Susy/MixingMatrix.h"
#include "ThePEG/PDT/EnumParticles.h"
#include "Herwig++/Looptools/clooptools.h"
using namespace ThePEG::Helicity;
using namespace Herwig;
namespace {
unsigned int neutralinoIndex(long id) {
if(id> 1000000)
return id<1000025 ? id-1000022 : (id-1000005)/10;
else if(abs(id)<=16)
return (abs(id)-4)/2;
else
return id-13;
}
}
SSGNGVertex::SSGNGVertex() : _includeOnShell(false), _realIntegral(false),
_omitLightQuarkYukawas(false),
_sw(0.), _cw(0.), _idlast(0),
_q2last(ZERO), _couplast(0.),
- _leftlast(ZERO), _rightlast(ZERO) {
+ _leftlast(ZERO), _rightlast(ZERO),
+ _initLoops(false) {
orderInGem(1);
orderInGs(2);
}
void SSGNGVertex::doinit() {
- Looptools::ltini();
+ if(!_initLoops) {
+ Looptools::ltini();
+ _initLoops = true;
+ }
tMSSMPtr model = dynamic_ptr_cast<tMSSMPtr>(generator()->standardModel());
if(!model)
throw InitException() << "SSGNGVertex::doinit() - "
<< "The model pointer is null."
<< Exception::abortnow;
_theN = model->neutralinoMix();
if(!_theN )
throw InitException() << "SSGNGVertex::doinit - The neutralino "
<< "mixing matrix pointer is null."
<< Exception::abortnow;
vector<long> ineu(4);
ineu[0] = 1000022; ineu[1] = 1000023;
ineu[2] = 1000025; ineu[3] = 1000035;
if(_theN->size().first==5)
ineu.push_back(1000045);
else if(_theN->size().first==7) {
if(model->majoranaNeutrinos()) {
ineu.push_back(17);
ineu.push_back(18);
ineu.push_back(19);
}
else {
ineu.push_back(12);
ineu.push_back(14);
ineu.push_back(16);
}
}
for(unsigned int i = 0; i < ineu.size(); ++i) {
addToList(1000021, ineu[i], 21);
}
GeneralFFVVertex::doinit();
_sw = sqrt(sin2ThetaW());
_cw = sqrt(1 - _sw*_sw);
_mw = getParticleData(ParticleID::Wplus)->mass();
double tb = model->tanBeta();
_sb = tb/sqrt(1 + sqr(tb));
_cb = sqrt(1 - sqr(_sb));
_stop = model->stopMix();
_sbot = model->sbottomMix();
Looptools::ltexi();
}
void SSGNGVertex::dofinish() {
Looptools::ltexi();
GeneralFFVVertex::dofinish();
}
void SSGNGVertex::doinitrun() {
- Looptools::ltini();
+ if(!_initLoops) {
+ Looptools::ltini();
+ _initLoops = true;
+ }
GeneralFFVVertex::doinitrun();
}
void SSGNGVertex::persistentOutput(PersistentOStream & os) const {
os << _sw << _cw << _theN << ounit(_mw,GeV) << _sb << _cb
<< _stop << _sbot << _includeOnShell << _omitLightQuarkYukawas
<< _realIntegral;
}
void SSGNGVertex::persistentInput(PersistentIStream & is, int) {
is >> _sw >> _cw >> _theN >> iunit(_mw,GeV) >> _sb >> _cb
>> _stop >> _sbot >> _includeOnShell >> _omitLightQuarkYukawas
>> _realIntegral;
}
// Static variable needed for the type description system in ThePEG.
DescribeClass<SSGNGVertex,GeneralFFVVertex>
describeHerwigSSGNGVertex("Herwig::SSGNGVertex", "HwSusy.so");
void SSGNGVertex::Init() {
static ClassDocumentation<SSGNGVertex> documentation
("The loop-mediated coupling of the gluino to a gluon and a neutralino.");
static Switch<SSGNGVertex,bool> interfaceIncludeOnShellIntermediates
("IncludeOnShellIntermediates",
"Whether or not to include on-shell intermediate states",
&SSGNGVertex::_includeOnShell, false, false, false);
static SwitchOption interfaceIncludeOnShellIntermediatesYes
(interfaceIncludeOnShellIntermediates,
"Yes",
"Include them",
true);
static SwitchOption interfaceIncludeOnShellIntermediatesNo
(interfaceIncludeOnShellIntermediates,
"No",
"Don't incldue them",
false);
static Switch<SSGNGVertex,bool> interfaceOmitLightQuarkYukawas
("OmitLightQuarkYukawas",
"Omit the yukawa type couplings for down, up, strange"
" and charm quarks, mainly for testing vs ISAJET",
&SSGNGVertex::_omitLightQuarkYukawas, false, false, false);
static SwitchOption interfaceOmitLightQuarkYukawasNo
(interfaceOmitLightQuarkYukawas,
"No",
"Include the Yukawas",
false);
static SwitchOption interfaceOmitLightQuarkYukawasYes
(interfaceOmitLightQuarkYukawas,
"Yes",
"Omit Yukawas",
true);
static Switch<SSGNGVertex,bool> interfaceRealIntegral
("RealIntegral",
"Only include the real parts of the integrals",
&SSGNGVertex::_realIntegral, false, false, false);
static SwitchOption interfaceRealIntegralYes
(interfaceRealIntegral,
"Yes",
"Only include the real part",
true);
static SwitchOption interfaceRealIntegralNo
(interfaceRealIntegral,
"No",
"Don't include the real part",
false);
}
void SSGNGVertex::setCoupling(Energy2 q2, tcPDPtr part1,
#ifndef NDEBUG
tcPDPtr part2,tcPDPtr part3) {
#else
tcPDPtr part2,tcPDPtr) {
#endif
int o[2]={1,0};
long in1 = part1->id();
long in2 = part2->id();
Energy Mj = part1->mass();
Energy Mi = part2->mass();
if(in1!=ParticleID::SUSY_g) {
swap(in1,in2);
swap(Mj,Mi);
}
// checks of the particle ids
assert(part3->id()==ParticleID::g);
assert(in1 == ParticleID::SUSY_g);
assert(in2 == ParticleID::SUSY_chi_10 || in2 == ParticleID::SUSY_chi_20 ||
in2 == ParticleID::SUSY_chi_30 || in2 == ParticleID::SUSY_chi_40 ||
in2 == 1000045 || in2==12 || in2==14 || in2==16 || in2==17 || in2==18 || in2==19 );
// normal couplings are zero
setLeft (0.);
setRight(0.);
if(in2 != _idlast || q2 !=_q2last) {
+ if(!_initLoops) {
+ Looptools::ltini();
+ _initLoops = true;
+ }
Looptools::clearcache();
_leftlast = ZERO;
_rightlast = ZERO;
_idlast = in2;
unsigned int neu = neutralinoIndex(in2);
Complex n1prime = (*_theN)(neu,0)*_cw + (*_theN)(neu,1)*_sw;
Complex n2prime = (*_theN)(neu,1)*_cw - (*_theN)(neu,0)*_sw;
// squark/quark loops
for(long iferm=1;iferm<7;++iferm) {
tcPDPtr smf = getParticleData(iferm);
Energy mf = smf->mass();
double qf = smf->charge()/eplus;
double y = (!(iferm<=4&&_omitLightQuarkYukawas)) ?
0.5*double(dynamic_ptr_cast<tcHwSMPtr>(generator()->standardModel())->mass(q2,smf)/_mw) : 0.;
Complex bracketl = qf*_sw*( conj(n1prime) - _sw*conj(n2prime)/_cw );
double lambda(0.);
// neutralino mixing element
Complex nlf(0.);
if( iferm % 2 == 0 ) {
y /= _sb;
lambda = -0.5 + qf*sqr(_sw);
nlf = (*_theN)(neu,3);
}
else {
y /= _cb;
lambda = 0.5 + qf*sqr(_sw);
nlf = (*_theN)(neu,2);
}
Complex bracketr = _sw*qf*n1prime - n2prime*lambda/_cw;
for(long iy=0;iy<2;++iy) {
long isf = 1000000*(1+iy)+iferm;
Energy msf = getParticleData(isf)->mass();
if(!_includeOnShell&&(mf+msf<Mj||mf+msf<Mi)) continue;
Complex g[2][2];
Complex ma1(0.), ma2(0.);
// heavy fermions
if( iferm == 5 || iferm == 6 ) {
if( iferm == 5 ) {
ma1 = (*_sbot)(iy,0);
ma2 = (*_sbot)(iy,1);
}
else if( iferm == 6 ) {
ma1 = (*_stop)(iy,0);
ma2 = (*_stop)(iy,1);
}
}
else if(iy==0) {
ma1 = 1.;
}
else {
ma2 = 1.;
}
g[0][0] = y*conj(nlf)*ma1 - ma2*bracketl;
g[0][1] = y*nlf *ma2 + ma1*bracketr;
g[1][0] = - ma2;
g[1][1] = + ma1;
swap(g[0][0],g[0][1]);
complex<InvEnergy2> I,J,K,I2;
loopIntegrals(Mi,Mj,msf,mf,I,J,K,I2);
complex<InvEnergy> coup[2];
for(unsigned int ix=0;ix<2;++ix) {
coup[ix] = Mj*(I2-K)*(g[0][ix]*g[1][o[ix]]-conj(g[0][o[ix]]*g[1][ix]))
+Mi*K*(g[0][o[ix]]*g[1][ix]-conj(g[0][ix]*g[1][o[ix]]))
+mf*I*(g[0][ix]*g[1][ix]-conj(g[0][o[ix]]*g[1][o[ix]]));
}
_leftlast += 2.*coup[0];
_rightlast += 2.*coup[1];
}
}
}
if(q2 != _q2last || _couplast==0.) {
_q2last = q2;
_couplast = weakCoupling(q2)*sqr(strongCoupling(q2))/
32./sqr(Constants::pi);
}
norm(_couplast);
setLeftSigma ( _leftlast);
- setRightSigma(_rightlast);
+ setRightSigma(_rightlast);
}
void SSGNGVertex::loopIntegrals(Energy Mi, Energy Mj, Energy M, Energy m,
complex<InvEnergy2> & I, complex<InvEnergy2> & J,
complex<InvEnergy2> & K, complex<InvEnergy2> & I2) {
Energy2 m2(sqr(m)),M2(sqr(M)),Mi2(sqr(Mi)),Mj2(sqr(Mj));
double min2 = Mj2*UnitRemoval::InvE2;
double mout2 = Mi2*UnitRemoval::InvE2;
double mf2 = m2 *UnitRemoval::InvE2;
double ms2 = M2 *UnitRemoval::InvE2;
I = Looptools::C0i(Looptools::cc0,min2,mout2,0.,mf2,ms2,mf2)*UnitRemoval::InvE2;
J = Looptools::C0i(Looptools::cc0,min2,mout2,0.,ms2,mf2,ms2)*UnitRemoval::InvE2;
I2 =-Looptools::C0i(Looptools::cc1,min2,mout2,0.,mf2,ms2,mf2)*UnitRemoval::InvE2;
K = (1.+Complex(m2*I+M2*J-Mj2*I2))/(Mi2-Mj2);
if(_realIntegral) {
I = I .real();
J = J .real();
I2 = I2.real();
K = K .real();
}
}
diff --git a/Models/Susy/SSGNGVertex.h b/Models/Susy/SSGNGVertex.h
--- a/Models/Susy/SSGNGVertex.h
+++ b/Models/Susy/SSGNGVertex.h
@@ -1,218 +1,223 @@
// -*- C++ -*-
//
// SSGNGVertex.h is a part of Herwig++ - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2011 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_SSGNGVertex_H
#define HERWIG_SSGNGVertex_H
//
// This is the declaration of the SSGNGVertex class.
//
#include "ThePEG/Helicity/Vertex/Vector/GeneralFFVVertex.h"
#include "Herwig++/Models/Susy/MSSM.h"
#include "Herwig++/Models/Susy/MixingMatrix.fh"
namespace Herwig {
using namespace ThePEG;
/**
* This class implements the coupling of a gluon to a gluino and a neutralino
* via loop diagrams
* It inherits from GeneralFFVertex and implements the setCoupling method.
*
* @see \ref SSGNGVertexInterfaces "The interfaces"
* defined for SSGNGVertex.
* @see FFVertex
*/
class SSGNGVertex: public GeneralFFVVertex {
public:
/**
* The default constructor.
*/
SSGNGVertex();
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();
/**
* Calculate the couplings.
* @param q2 The scale \f$q^2\f$ for the coupling at the vertex.
* @param part1 The ParticleData pointer for the first particle.
* @param part2 The ParticleData pointer for the second particle.
* @param part3 The ParticleData pointer for the third particle.
*/
virtual void setCoupling(Energy2 q2, tcPDPtr part1,
tcPDPtr part2, tcPDPtr part3);
protected:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const {return new_ptr(*this);}
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const {return new_ptr(*this);}
//@}
protected:
/**
* Evaluate the loop integrals
*/
void loopIntegrals(Energy Mi, Energy Mj, Energy M, Energy m,
complex<InvEnergy2> & I, complex<InvEnergy2> & J,
complex<InvEnergy2> & K, complex<InvEnergy2> & I2);
protected:
/** @name Standard Interfaced functions. */
//@{
/**
* Initialize this object after the setup phase before saving an
* EventGenerator to disk.
* @throws InitException if object could not be initialized properly.
*/
virtual void doinit();
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
/**
* Finalize this object. Called in the run phase just after a
* run has ended. Used eg. to write out statistics.
*/
virtual void dofinish();
//@}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
SSGNGVertex & operator=(const SSGNGVertex &);
private:
/**
* Whether of not to include on-shell intermediate states
*/
bool _includeOnShell;
/**
* Only include the real part of the integral
*/
bool _realIntegral;
/**
* Option to omit light quark yukawas
*/
bool _omitLightQuarkYukawas;
/**
* Pointer to the stop mixing matrix
*/
tMixingMatrixPtr _stop;
/**
* Pointer to the sbottom mixing matrix
*/
tMixingMatrixPtr _sbot;
/**
* The value of \f$sin(\theta_W)\f$
*/
double _sw;
/**
* The value of \f$cos(\theta_W)\f$
*/
double _cw;
/**
* Mass of the W
*/
Energy _mw;
/**
* \f$\sin(\beta)\f$
*/
double _sb;
/**
* \f$\cos(\beta)\f$
*/
double _cb;
/**
* Store the neutralino mixing matrix
*/
tMixingMatrixPtr _theN;
/**
* Store the id of the neutralino when the coupling was last evaluated
*/
long _idlast;
/**
- * Store the value at which the coupling when it was last evalutated
+ * Store the value at which the coupling when it was last evaluated
*/
Energy2 _q2last;
/**
- * Store the value of the coupling when it was last evalutated
+ * Store the value of the coupling when it was last evaluated
*/
Complex _couplast;
/**
- * Store the value of the left-coupling when it was last evalutated
+ * Store the value of the left-coupling when it was last evaluated
*/
complex<InvEnergy> _leftlast;
/**
- * Store the value of the right-coupling when it was last evalutated
+ * Store the value of the right-coupling when it was last evaluated
*/
complex<InvEnergy> _rightlast;
+
+ /**
+ * Whether or not initialised
+ */
+ bool _initLoops;
};
}
#endif /* HERWIG_SSGNGVertex_H */
diff --git a/Models/TTbAsymm/TTbAModel.h b/Models/TTbAsymm/TTbAModel.h
--- a/Models/TTbAsymm/TTbAModel.h
+++ b/Models/TTbAsymm/TTbAModel.h
@@ -1,354 +1,354 @@
// -*- C++ -*-
#ifndef HERWIG_TTbAModel_H
#define HERWIG_TTbAModel_H
//
// This is the declaration of the TTbAModel class.
//
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "TTbAModel.fh"
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
/**
* Here is the documentation of the TTbAModel class.
*
* @see \ref TTbAModelInterfaces "The interfaces"
* defined for TTbAModel.
*/
-class TTbAModel: public StandardModel {
+class TTbAModel: public BSMModel {
public:
/**
* The default constructor.
*/
TTbAModel();
/** @name Vertices */
//@{
/**
* Pointer to the object handling W prime vertex.
*/
tAbstractFFVVertexPtr vertexWPTD() const {return _theWPTDVertex;}
/**
* Pointer to the object handling Z prime vertex.
*/
tAbstractFFVVertexPtr vertexZPQQ() const {return _theZPQQVertex;}
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();
/**
* Return the W prime top-down left-handed coupling
*/
double _cWPTD_left() const {return _gWPTD_L;}
/**
* Return the W prime top-down right-handed coupling
*/
double _cWPTD_right() const {return _gWPTD_R;}
/**
* Return the Z prime top-up left-handed coupling
*/
double _cZPTU_left() const {return _gZPTU_L;}
/**
* Return the Z prime top-up right-handed coupling
*/
double _cZPTU_right() const {return _gZPTU_R;}
/**
* Return the Z prime up-upbar left-handed coupling
*/
double _cZPUU_left() const {return _gZPUU_L;}
/**
* Return the Z prime up-upbar right-handed coupling
*/
double _cZPUU_right() const {return _gZPUU_R;}
/**
* Return the Z prime charm-charmbar left-handed coupling
*/
double _cZPCC_left() const {return _gZPCC_L;}
/**
* Return the Z prime charm-charmbar right-handed coupling
*/
double _cZPCC_right() const {return _gZPCC_R;}
/**
* Return the axigluon q-qbar left-handed coupling
*/
double _cAGQQ_left() const {return _gAGQQ_L;}
/**
* Return the axigluon q-qbar right-handed coupling
*/
double _cAGQQ_right() const {return _gAGQQ_R;}
/**
* Return the axigluon t-tbar left-handed coupling
*/
double _cAGTT_left() const {return _gAGTT_L;}
/**
* Return the axigluon t-tbar right-handed coupling
*/
double _cAGTT_right() const {return _gAGTT_R;}
/**
* Return the alphaX value of the SU(2)_X model
*/
double _alphaX_value() const {return _alphaXparam;}
/**
* Return the costheta misalignment value of the SU(2)_X model
*/
double _costhetaX_value() const {return _costhetaXparam;}
/**
* Return the selected model id
*/
int _model() const {return _modelselect;}
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:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<TTbAModel> initTTbAModel;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
TTbAModel & operator=(const TTbAModel &);
/**
* Pointer to the object handling the Wp to Top Down vertex.
*/
AbstractFFVVertexPtr _theWPTDVertex;
/**
* Pointer to the object handling the Zp to Quark-antiQuark vertex.
*/
AbstractFFVVertexPtr _theZPQQVertex;
/**
* Pointer to the object handling the Ag to Quark-antiQuark vertex.
*/
AbstractFFVVertexPtr _theAGQQVertex;
/**
* Pointer to the object handling the SU(2)_X vertex.
*/
AbstractFFVVertexPtr _theSU2XVertex;
/**
* W prime coupling to top-down (left-handed)
*/
double _gWPTD_L;
/**
* W prime coupling to top-down (right-handed)
*/
double _gWPTD_R;
/**
* Z prime coupling to top-up (left-handed)
*/
double _gZPTU_L;
/**
* Z prime coupling to top-up (right-handed)
*/
double _gZPTU_R;
/**
* Z prime coupling to up-upbar (left-handed)
*/
double _gZPUU_L;
/**
* Z prime coupling to up-upbar (right-handed)
*/
double _gZPUU_R;
/**
* Z prime coupling to charm-charmbar (left-handed)
*/
double _gZPCC_L;
/**
* Z prime coupling to charm-charmbar (right-handed)
*/
double _gZPCC_R;
/**
* Axigluon coupling to q-qbar (left-handed)
*/
double _gAGQQ_L;
/**
* Axigluon coupling to q-qbar (right-handed)
*/
double _gAGQQ_R;
/**
* Axigluon coupling to t-tbar (left-handed)
*/
double _gAGTT_L;
/**
* Axigluon coupling to t-tbar (right-handed)
*/
double _gAGTT_R;
/**
* SU(2)_X alpha_X parameter
*/
double _alphaXparam;
/**
* SU(2)_X costheta misalignment angle
*/
double _costhetaXparam;
/**
* Model selector
*/
int _modelselect;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of TTbAModel. */
template <>
struct BaseClassTrait<Herwig::TTbAModel,1> {
/** Typedef of the first base class of TTbAModel. */
typedef Herwig::StandardModel NthBase;
};
/** This template specialization informs ThePEG about the name of
* the TTbAModel class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::TTbAModel>
: public ClassTraitsBase<Herwig::TTbAModel> {
/** Return a platform-independent class name */
static string className() { return "Herwig::TTbAModel"; }
/**
* The name of a file containing the dynamic library where the class
* TTbAModel is implemented. It may also include several, space-separated,
* libraries if the class TTbAModel depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwTTbAModel.so"; }
};
/** @endcond */
}
#endif /* HERWIG_TTbAModel_H */
diff --git a/Models/Zprime/ZprimeModel.h b/Models/Zprime/ZprimeModel.h
--- a/Models/Zprime/ZprimeModel.h
+++ b/Models/Zprime/ZprimeModel.h
@@ -1,465 +1,465 @@
// -*- C++ -*-
#ifndef HERWIG_ZprimeModel_H
#define HERWIG_ZprimeModel_H
//
// This is the declaration of the ZprimeModel class.
//
-#include "Herwig++/Models/StandardModel/StandardModel.h"
+#include "Herwig++/Models/General/BSMModel.h"
#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.h"
#include "ZprimeModel.fh"
namespace Herwig {
using namespace ThePEG;
using namespace ThePEG::Helicity;
/**
* Here is the documentation of the ZprimeModel class.
*
* @see \ref ZprimeModelInterfaces "The interfaces"
* defined for ZprimeModel.
*/
-class ZprimeModel: public StandardModel {
+class ZprimeModel: public BSMModel {
public:
/**
* The default constructor.
*/
ZprimeModel();
/** @name Vertices */
//@{
/**
* Pointer to the object handling Z prime quark-anti-quark vertex.
*/
tAbstractFFVVertexPtr vertexZPQQ() const {return _theZPQQVertex;}
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();
/**
* Return the Z prime top-up left-handed coupling
*/
double _cZPTU_left() const {return _gZPTU_L;}
/**
* Return the Z prime top-up right-handed coupling
*/
double _cZPTU_right() const {return _gZPTU_R;}
/**
* Return the Z prime d-dbar left-handed coupling
*/
double _cZPDD_left() const {return _gZPDD_L;}
/**
* Return the Z prime d-dbar right-handed coupling
*/
double _cZPDD_right() const {return _gZPDD_R;}
/**
* Return the Z prime top-anti-top left-handed coupling
*/
double _cZPTT_left() const {return _gZPTT_L;}
/**
* Return the Z prime top-anti-top right-handed coupling
*/
double _cZPTT_right() const {return _gZPTT_R;}
/**
* Return the Z prime u-ubar left-handed coupling
*/
double _cZPUU_left() const {return _gZPUU_L;}
/**
* Return the Z prime u-ubar right-handed coupling
*/
double _cZPUU_right() const {return _gZPUU_R;}
/**
* Return the Z prime c-cbar left-handed coupling
*/
double _cZPCC_left() const {return _gZPCC_L;}
/**
* Return the Z prime c-cbar right-handed coupling
*/
double _cZPCC_right() const {return _gZPCC_R;}
/**
* Return the Z prime b-bbar left-handed coupling
*/
double _cZPBB_left() const {return _gZPBB_L;}
/**
* Return the Z prime b-bbar right-handed coupling
*/
double _cZPBB_right() const {return _gZPBB_R;}
/**
* Return the Z prime s-sbar left-handed coupling
*/
double _cZPSS_left() const {return _gZPSS_L;}
/**
* Return the Z prime c-cbar right-handed coupling
*/
double _cZPSS_right() const {return _gZPSS_R;}
/**
* Return the Z prime e+e- left-handed coupling
*/
double _cZPee_left() const {return _gZPee_L;}
/**
* Return the Z prime e+e- right-handed coupling
*/
double _cZPee_right() const {return _gZPee_R;}
/**
* Return the Z prime mu+mu- left-handed coupling
*/
double _cZPmm_left() const {return _gZPmm_L;}
/**
* Return the Z prime mu+mu- right-handed coupling
*/
double _cZPmm_right() const {return _gZPmm_R;}
/**
* Return the Z prime tau+tau- left-handed coupling
*/
double _cZPtt_left() const {return _gZPtt_L;}
/**
* Return the Z prime tau+tau- right-handed coupling
*/
double _cZPtt_right() const {return _gZPtt_R;}
/**
* Return the Z prime nu_e nu_ebar left-handed coupling
*/
double _cZPnuenue_left() const {return _gZPnuenue_L;}
/**
* Return the Z prime nu_e nu_ebar right-handed coupling
*/
double _cZPnuenue_right() const {return _gZPnuenue_R;}
/**
* Return the Z prime nu_mu nu_mubar left-handed coupling
*/
double _cZPnumnum_left() const {return _gZPnumnum_L;}
/**
* Return the Z prime nu_mu nu_mubar right-handed coupling
*/
double _cZPnumnum_right() const {return _gZPnumnum_R;}
/**
* Return the Z prime nu_tau nu_taubar left-handed coupling
*/
double _cZPnutnut_left() const {return _gZPnutnut_L;}
/**
* Return the Z prime nu_tau nu_taubar right-handed coupling
*/
double _cZPnutnut_right() const {return _gZPnutnut_R;}
/**
* Return the overall coupling of the Z prime to quark-anti-quark
*/
double _cZPoverallCoup() const {return _ZPoverall;}
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:
/** @name Clone Methods. */
//@{
/**
* Make a simple clone of this object.
* @return a pointer to the new object.
*/
virtual IBPtr clone() const;
/** Make a clone of this object, possibly modifying the cloned object
* to make it sane.
* @return a pointer to the new object.
*/
virtual IBPtr fullclone() const;
//@}
// If needed, insert declarations of virtual function defined in the
// InterfacedBase class here (using ThePEG-interfaced-decl in Emacs).
private:
/**
* The static object used to initialize the description of this class.
* Indicates that this is a concrete class with persistent data.
*/
static ClassDescription<ZprimeModel> initZprimeModel;
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
ZprimeModel & operator=(const ZprimeModel &);
/**
* Pointer to the object handling the Zp to Quark-antiQuark vertex.
*/
AbstractFFVVertexPtr _theZPQQVertex;
/**
* Z prime coupling to u-ubar (left-handed)
*/
double _gZPUU_L;
/**
* Z prime coupling to u-ubar (right-handed)
*/
double _gZPUU_R;
/**
* Z prime coupling to d-dbar (left-handed)
*/
double _gZPDD_L;
/**
* Z prime coupling to d-dbar (right-handed)
*/
double _gZPDD_R;
/**
* Z prime coupling to c-cbar (left-handed)
*/
double _gZPCC_L;
/**
* Z prime coupling to c-cbar (right-handed)
*/
double _gZPCC_R;
/**
* Z prime coupling to s-sbar (left-handed)
*/
double _gZPSS_L;
/**
* Z prime coupling to s-sbar (right-handed)
*/
double _gZPSS_R;
/**
* Z prime coupling to b-bbar (left-handed)
*/
double _gZPBB_L;
/**
* Z prime coupling to b-bbar (right-handed)
*/
double _gZPBB_R;
/**
* Z prime coupling to top-up (left-handed)
*/
double _gZPTU_L;
/**
* Z prime coupling to top-up (right-handed)
*/
double _gZPTU_R;
/**
* Z prime coupling to top-anti-top (left-handed)
*/
double _gZPTT_L;
/**
* Z prime coupling to top-anti-top (right-handed)
*/
double _gZPTT_R;
/**
* Z prime coupling to e+e- (left-handed)
*/
double _gZPee_L;
/**
* Z prime coupling to e+e- (right-handed)
*/
double _gZPee_R;
/**
* Z prime coupling to mu+mu- (left-handed)
*/
double _gZPmm_L;
/**
* Z prime coupling to mu+mu- (right-handed)
*/
double _gZPmm_R;
/**
* Z prime coupling to tau+tau- (left-handed)
*/
double _gZPtt_L;
/**
* Z prime coupling to tau+tau- (right-handed)
*/
double _gZPtt_R;
/**
* Z prime coupling to nu_e nu_ebar (left-handed)
*/
double _gZPnuenue_L;
/**
* Z prime coupling to nu_e nu_ebar (right-handed)
*/
double _gZPnuenue_R;
/**
* Z prime coupling to nu_mu nu_mubar (left-handed)
*/
double _gZPnumnum_L;
/**
* Z prime coupling to nu_mu nu_mubar (right-handed)
*/
double _gZPnumnum_R;
/**
* Z prime coupling to nu_tau nu_taubar (left-handed)
*/
double _gZPnutnut_L;
/**
* Z prime coupling to nu_tau nu_taubar (right-handed)
*/
double _gZPnutnut_R;
/**
* Z prime overall coupling
*/
double _ZPoverall;
};
}
#include "ThePEG/Utilities/ClassTraits.h"
namespace ThePEG {
/** @cond TRAITSPECIALIZATIONS */
/** This template specialization informs ThePEG about the
* base classes of ZprimeModel. */
template <>
struct BaseClassTrait<Herwig::ZprimeModel,1> {
/** Typedef of the first base class of ZprimeModel. */
typedef Herwig::StandardModel NthBase;
};
/** This template specialization informs ThePEG about the name of
* the ZprimeModel class and the shared object where it is defined. */
template <>
struct ClassTraits<Herwig::ZprimeModel>
: public ClassTraitsBase<Herwig::ZprimeModel> {
/** Return a platform-independent class name */
static string className() { return "Herwig::ZprimeModel"; }
/**
* The name of a file containing the dynamic library where the class
* ZprimeModel is implemented. It may also include several, space-separated,
* libraries if the class ZprimeModel depends on other classes (base classes
* excepted). In this case the listed libraries will be dynamically
* linked in the order they are specified.
*/
static string library() { return "HwZprimeModel.so"; }
};
/** @endcond */
}
#endif /* HERWIG_ZprimeModel_H */
diff --git a/src/LHC-MB.in b/src/LHC-MB.in
--- a/src/LHC-MB.in
+++ b/src/LHC-MB.in
@@ -1,58 +1,54 @@
################################################################################
# This file contains our best tune to UE data from ATLAS at 7 TeV. More recent
# tunes and tunes for other centre-of-mass energies as well as more usage
# instructions can be obtained from this Herwig++ wiki page:
# http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes
################################################################################
##################################################
# Technical parameters for this run
##################################################
cd /Herwig/Generators
set LHCGenerator:NumberOfEvents 10000000
set LHCGenerator:RandomNumberGenerator:Seed 31122001
set LHCGenerator:PrintEvent 10
set LHCGenerator:MaxErrors 1000000
set LHCGenerator:DebugLevel 0
set LHCGenerator:DumpPeriod -1
set LHCGenerator:DebugEvent 0
##################################################
# LHC physics parameters (override defaults here)
##################################################
set LHCGenerator:EventHandler:LuminosityFunction:Energy 8000.0
# Intrinsic pT tune extrapolated to LHC energy
set /Herwig/Shower/Evolver:IntrinsicPtGaussian 2.2*GeV
##################################################
# Matrix Elements for hadron-hadron collisions
##################################################
cd /Herwig/MatrixElements/
insert SimpleQCD:MatrixElements[0] MEMinBias
-# Need this cut only for min bias
-cd /Herwig/Cuts
-set JetKtCut:MinKT 0.0*GeV
-set QCDCuts:MHatMin 0.0*GeV
-set QCDCuts:X1Min 0.055
-set QCDCuts:X2Min 0.055
# MPI model settings
set /Herwig/UnderlyingEvent/MPIHandler:IdenticalToUE 0
+cd /Herwig/Generators
-cd /Herwig/Generators
+# Change to have no pT cuts for MinBias
+set LHCGenerator:EventHandler:Cuts /Herwig/Cuts/MinBiasCuts
#insert LHCGenerator:AnalysisHandlers 0 /Herwig/Analysis/HepMCFile
#set /Herwig/Analysis/HepMCFile:PrintEvent 1000000
#set /Herwig/Analysis/HepMCFile:Format GenEvent
#set /Herwig/Analysis/HepMCFile:Units GeV_mm
#set /Herwig/Analysis/HepMCFile:Filename events.fifo
##################################################
# Save run for later usage with 'Herwig++ run'
##################################################
saverun LHC-MB LHCGenerator
diff --git a/src/defaults/BSM.in b/src/defaults/BSM.in
--- a/src/defaults/BSM.in
+++ b/src/defaults/BSM.in
@@ -1,126 +1,128 @@
#################################################
# Create the BSM machinery
# This only gets switched on in read if it is
# requested by the user
################################################
mkdir /Herwig/NewPhysics
cd /Herwig/NewPhysics
create Herwig::ModelGenerator NewModel
create Herwig::ResonantProcessConstructor ResConstructor
create Herwig::TwoToTwoProcessConstructor HPConstructor
create Herwig::HiggsVectorBosonProcessConstructor HVConstructor
set HVConstructor:AlphaQCD /Herwig/Shower/AlphaQCD
create Herwig::HiggsVBFProcessConstructor HiggsVBFConstructor
create Herwig::QQHiggsProcessConstructor QQHiggsConstructor
create Herwig::DecayConstructor DecayConstructor
newdef DecayConstructor:QEDGenerator /Herwig/QEDRadiation/SOPHTY
create Herwig::TwoBodyDecayConstructor TwoBodyDC
set TwoBodyDC:ExcludeEffectiveVertices No
create Herwig::ThreeBodyDecayConstructor ThreeBodyDC
+set ThreeBodyDC:RemoveFlavourChangingVertices Yes
create Herwig::FourBodyDecayConstructor FourBodyDC
+set FourBodyDC:RemoveFlavourChangingVertices Yes
create Herwig::WeakCurrentDecayConstructor WeakDecayConstructor
set WeakDecayConstructor:InitPoints 10000
set WeakDecayConstructor:MassCut 2.
# pi-
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau1MesonCurrent
insert WeakDecayConstructor:DecayModes 0 pi-;
insert WeakDecayConstructor:Normalisation 0 1.01386262897
# pi-,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2MesonCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi0;
insert WeakDecayConstructor:Normalisation 0 1.17616809738
# e-,nu_ebar
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2LeptonCurrent
insert WeakDecayConstructor:DecayModes 0 e-,nu_ebar;
insert WeakDecayConstructor:Normalisation 0 1.
# mu-,nu_mubar
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2LeptonCurrent
insert WeakDecayConstructor:DecayModes 0 mu-,nu_mubar;
insert WeakDecayConstructor:Normalisation 0 1.
# tau-,nu_taubar
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2LeptonCurrent
insert WeakDecayConstructor:DecayModes 0 tau-,nu_taubar;
insert WeakDecayConstructor:Normalisation 0 1.
# pi-,pi0,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi0,pi0;
insert WeakDecayConstructor:Normalisation 0 1.65956712121
# pi-,pi+,pi-
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi+,pi-;
insert WeakDecayConstructor:Normalisation 0 1.62175791702
# pi-,pi+,pi-,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau4PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi+,pi-,pi0;
insert WeakDecayConstructor:Normalisation 0 1.09170097618
# pi-,pi0,pi0,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau4PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi0,pi0,pi0;
insert WeakDecayConstructor:Normalisation 0 1.
# Kbar0,pi-
insert WeakDecayConstructor:Current 0 /Herwig/Decays/TauKPiCurrent
insert WeakDecayConstructor:DecayModes 0 Kbar0,pi-;
insert WeakDecayConstructor:Normalisation 0 1.17156865176
# K-
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau1MesonCurrent
insert WeakDecayConstructor:DecayModes 0 K-;
insert WeakDecayConstructor:Normalisation 0 0.938794564668
# K-,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/TauKPiCurrent
insert WeakDecayConstructor:DecayModes 0 K-,pi0;
insert WeakDecayConstructor:Normalisation 0 1.12526014943
# pi-,pi-,pi+,pi0,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau5PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi-,pi-,pi+,pi0,pi0;
insert WeakDecayConstructor:Normalisation 0 0.954286993254
# pi-,pi-,pi-,pi+,pi+
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau5PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi+,pi+,pi+,pi-,pi-;
insert WeakDecayConstructor:Normalisation 0 0.435953860245
# pi-,pi0,pi0,pi0,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau5PionCurrent
insert WeakDecayConstructor:DecayModes 0 pi+,pi0,pi0,pi0,pi0;
insert WeakDecayConstructor:Normalisation 0 0.603378959531
# K0,pi+,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3KaonCurrent
insert WeakDecayConstructor:DecayModes 0 K0,pi+,pi0;
insert WeakDecayConstructor:Normalisation 0 0.380268556539
# K+,pi+,pi-
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3KaonCurrent
insert WeakDecayConstructor:DecayModes 0 K+,pi+,pi-;
insert WeakDecayConstructor:Normalisation 0 0.407904176498
# K+,K-,pi+
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3KaonCurrent
insert WeakDecayConstructor:DecayModes 0 K+,K-,pi+;
insert WeakDecayConstructor:Normalisation 0 0.727416124384
# K+,Kbar0,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3KaonCurrent
insert WeakDecayConstructor:DecayModes 0 K+,Kbar0,pi0;
insert WeakDecayConstructor:Normalisation 0 0.945076580717
# K_L0,K_S0,pi+
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3KaonCurrent
insert WeakDecayConstructor:DecayModes 0 K_L0,K_S0,pi+;
insert WeakDecayConstructor:Normalisation 0 1.10729951668
# pi+,pi0,gamma
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2MesonPhotonCurrent
insert WeakDecayConstructor:DecayModes 0 pi+,pi0,gamma;
insert WeakDecayConstructor:Normalisation 0 3.13799210805
# eta,pi+,pi0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau3MesonCurrent
insert WeakDecayConstructor:DecayModes 0 eta,pi+,pi0;
insert WeakDecayConstructor:Normalisation 0 0.880539616307
# K+,Kbar0
insert WeakDecayConstructor:Current 0 /Herwig/Decays/Tau2MesonCurrent
insert WeakDecayConstructor:DecayModes 0 K+,Kbar0;
insert WeakDecayConstructor:Normalisation 0 1.62501804004
insert NewModel:HardProcessConstructors[0] HPConstructor
insert NewModel:HardProcessConstructors[1] ResConstructor
insert NewModel:HardProcessConstructors[2] HVConstructor
insert NewModel:HardProcessConstructors[3] HiggsVBFConstructor
insert NewModel:HardProcessConstructors[4] QQHiggsConstructor
newdef NewModel:DecayConstructor DecayConstructor
insert DecayConstructor:NBodyDecayConstructors[0] TwoBodyDC
insert DecayConstructor:NBodyDecayConstructors[1] ThreeBodyDC
insert DecayConstructor:NBodyDecayConstructors[2] FourBodyDC
diff --git a/src/defaults/Cuts.in b/src/defaults/Cuts.in
--- a/src/defaults/Cuts.in
+++ b/src/defaults/Cuts.in
@@ -1,115 +1,122 @@
###########################################################
# Default cuts (applied to the hard subprocess)
#
# Don't change values here, re-set them in your own input
# files using these as examples.
###########################################################
mkdir /Herwig/Matchers
cd /Herwig/Matchers
create ThePEG::Matcher<Lepton> Lepton
create ThePEG::Matcher<ChargedLepton> ChargedLepton
create ThePEG::Matcher<LightQuark> LightQuark
create ThePEG::Matcher<LightAntiQuark> LightAntiQuark
create ThePEG::Matcher<StandardQCDParton> StandardQCDParton
create ThePEG::Matcher<Photon> Photon
create ThePEG::Matcher<Top> Top
create ThePEG::Matcher<WBoson> WBoson
create ThePEG::Matcher<ZBoson> ZBoson
create ThePEG::Matcher<HiggsBoson> HiggsBoson
mkdir /Herwig/Cuts
cd /Herwig/Cuts
# create the cuts object for e+e-
create ThePEG::Cuts EECuts
newdef EECuts:MHatMin 22.36*GeV
# create the cuts object for hadron collisions
create ThePEG::Cuts QCDCuts
newdef QCDCuts:ScaleMin 2.0*GeV
newdef QCDCuts:X1Min 1.0e-5
newdef QCDCuts:X2Min 1.0e-5
newdef QCDCuts:MHatMin 20.*GeV
# cut on jet pt
create ThePEG::SimpleKTCut JetKtCut SimpleKTCut.so
newdef JetKtCut:Matcher /Herwig/Matchers/StandardQCDParton
newdef JetKtCut:MinKT 20.0*GeV
# cut on photon
create ThePEG::SimpleKTCut PhotonKtCut SimpleKTCut.so
newdef PhotonKtCut:Matcher /Herwig/Matchers/Photon
newdef PhotonKtCut:MinKT 20.0*GeV
newdef PhotonKtCut:MinEta -3.
newdef PhotonKtCut:MaxEta 3.
# cut on leptons
create ThePEG::SimpleKTCut LeptonKtCut SimpleKTCut.so
newdef LeptonKtCut:Matcher /Herwig/Matchers/Lepton
newdef LeptonKtCut:MinKT 0.0*GeV
# cut on charged leptons
create ThePEG::SimpleKTCut ChargedLeptonKtCut SimpleKTCut.so
newdef ChargedLeptonKtCut:Matcher /Herwig/Matchers/ChargedLepton
newdef ChargedLeptonKtCut:MinKT 0.0*GeV
# cut on top quarks
create ThePEG::KTRapidityCut TopKtCut KTRapidityCut.so
newdef TopKtCut:Matcher /Herwig/Matchers/Top
newdef TopKtCut:MinKT 0.0*GeV
# cut on W bosons
create ThePEG::KTRapidityCut WBosonKtCut KTRapidityCut.so
newdef WBosonKtCut:Matcher /Herwig/Matchers/WBoson
newdef WBosonKtCut:MinKT 0.0*GeV
# cut on Z bosons
create ThePEG::KTRapidityCut ZBosonKtCut KTRapidityCut.so
newdef ZBosonKtCut:Matcher /Herwig/Matchers/ZBoson
newdef ZBosonKtCut:MinKT 0.0*GeV
# cut on Higgs bosons
create ThePEG::KTRapidityCut HiggsBosonKtCut KTRapidityCut.so
newdef HiggsBosonKtCut:Matcher /Herwig/Matchers/HiggsBoson
newdef HiggsBosonKtCut:MinKT 0.0*GeV
# create a cut on the invariant mass of lepton pairs
create ThePEG::V2LeptonsCut MassCut V2LeptonsCut.so
newdef MassCut:Families All
newdef MassCut:CComb All
newdef MassCut:MinM 20.*GeV
newdef MassCut:MaxM 14000.*GeV
# create a cut on Q^2 for neutral current DIS
create ThePEG::SimpleDISCut NeutralCurrentCut SimpleDISCut.so
newdef NeutralCurrentCut:MinQ2 20.
newdef NeutralCurrentCut:Current Neutral
# create a cut on Q^2 for charged current DIS
create ThePEG::SimpleDISCut ChargedCurrentCut SimpleDISCut.so
newdef ChargedCurrentCut:MinQ2 20.
newdef ChargedCurrentCut:Current Charged
# create a cut of Q^2 for charged current DIS
# insert into hadron cuts
insert QCDCuts:OneCuts[0] JetKtCut
insert QCDCuts:OneCuts[1] PhotonKtCut
insert QCDCuts:OneCuts[2] LeptonKtCut
insert QCDCuts:OneCuts[3] TopKtCut
insert QCDCuts:OneCuts[4] WBosonKtCut
insert QCDCuts:OneCuts[5] ZBosonKtCut
insert QCDCuts:OneCuts[6] HiggsBosonKtCut
insert QCDCuts:OneCuts[7] ChargedLeptonKtCut
insert QCDCuts:MultiCuts[0] MassCut
# cuts for DIS
create ThePEG::Cuts DISCuts
newdef DISCuts:ScaleMin 1.0*GeV
newdef DISCuts:X1Min 1.0e-5
newdef DISCuts:X2Min 1.0e-5
insert DISCuts:TwoCuts[0] NeutralCurrentCut
insert DISCuts:TwoCuts[1] ChargedCurrentCut
+
+# create diffrent cuts object for MinBias to avoid numerical problems
+create ThePEG::Cuts MinBiasCuts
+newdef MinBiasCuts:ScaleMin 2.0*GeV
+newdef MinBiasCuts:X1Min 0.055
+newdef MinBiasCuts:X2Min 0.055
+newdef MinBiasCuts:MHatMin 0.0*GeV

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 9:11 PM (1 d, 29 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3806213
Default Alt Text
(275 KB)

Event Timeline