diff --git a/Models/General/HardProcessConstructor.cc b/Models/General/HardProcessConstructor.cc
--- a/Models/General/HardProcessConstructor.cc
+++ b/Models/General/HardProcessConstructor.cc
@@ -1,835 +1,829 @@
// -*- C++ -*-
//
// This is the implementation of the non-inlined, non-templated member
// functions of the HardProcessConstructor class.
//
#include "HardProcessConstructor.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Interface/Switch.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
using namespace Herwig;
void HardProcessConstructor::persistentOutput(PersistentOStream & os) const {
os << debug_ << subProcess_ << model_;
}
void HardProcessConstructor::persistentInput(PersistentIStream & is, int) {
is >> debug_ >> subProcess_ >> model_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeAbstractClass<HardProcessConstructor,Interfaced>
describeHerwigHardProcessConstructor("Herwig::HardProcessConstructor", "Herwig.so");
void HardProcessConstructor::Init() {
static ClassDocumentation<HardProcessConstructor> documentation
("Base class for implementation of the automatic generation of hard processes");
static Switch<HardProcessConstructor,bool> interfaceDebugME
("DebugME",
"Print comparison with analytical ME",
&HardProcessConstructor::debug_, false, false, false);
static SwitchOption interfaceDebugMEYes
(interfaceDebugME,
"Yes",
"Print the debug information",
true);
static SwitchOption interfaceDebugMENo
(interfaceDebugME,
"No",
"Do not print the debug information",
false);
}
void HardProcessConstructor::doinit() {
Interfaced::doinit();
EGPtr eg = generator();
model_ = dynamic_ptr_cast<HwSMPtr>(eg->standardModel());
if(!model_)
throw InitException() << "HardProcessConstructor:: doinit() - "
<< "The model pointer is null!"
<< Exception::abortnow;
if(!eg->eventHandler()) {
throw
InitException() << "HardProcessConstructor:: doinit() - "
<< "The eventHandler pointer was null therefore "
<< "could not get SubProcessHandler pointer "
<< Exception::abortnow;
}
string subProcessName =
eg->preinitInterface(eg->eventHandler(), "SubProcessHandlers", "get","");
subProcess_ = eg->getObject<SubProcessHandler>(subProcessName);
if(!subProcess_) {
ostringstream s;
s << "HardProcessConstructor:: doinit() - "
<< "There was an error getting the SubProcessHandler "
<< "from the current event handler. ";
generator()->logWarning( Exception(s.str(), Exception::warning) );
}
}
GeneralHardME::ColourStructure HardProcessConstructor::
colourFlow(const tcPDVector & extpart) const {
PDT::Colour ina = extpart[0]->iColour();
PDT::Colour inb = extpart[1]->iColour();
PDT::Colour outa = extpart[2]->iColour();
PDT::Colour outb = extpart[3]->iColour();
// incoming colour neutral
if(ina == PDT::Colour0 && inb == PDT::Colour0) {
if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour11to11;
}
else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
return GeneralHardME::Colour11to33bar;
}
else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
return GeneralHardME::Colour11to88;
}
else
assert(false);
}
// incoming 3 3
else if(ina == PDT::Colour3 && inb == PDT::Colour3 ) {
if( outa == PDT::Colour3 && outb == PDT::Colour3 ) {
return GeneralHardME::Colour33to33;
}
else if( outa == PDT::Colour6 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour33to61;
}
else if( outa == PDT::Colour0 && outb == PDT::Colour6 ) {
return GeneralHardME::Colour33to16;
}
else if ( outa == PDT::Colour0 && outb == PDT::Colour3bar) {
return GeneralHardME::Colour33to13bar;
}
else if ( outb == PDT::Colour0 && outa == PDT::Colour3bar) {
return GeneralHardME::Colour33to3bar1;
}
else if ( outa == PDT::Colour8 && outb == PDT::Colour3bar) {
return GeneralHardME::Colour33to83bar;
}
else if ( outb == PDT::Colour8 && outa == PDT::Colour3bar) {
return GeneralHardME::Colour33to3bar8;
}
else
assert(false);
}
// incoming 3bar 3bar
else if(ina == PDT::Colour3bar && inb == PDT::Colour3bar ) {
if( outa == PDT::Colour3bar && outb == PDT::Colour3bar ) {
return GeneralHardME::Colour3bar3barto3bar3bar;
}
else if( outa == PDT::Colour6bar && outb == PDT::Colour0) {
return GeneralHardME::Colour3bar3barto6bar1;
}
else if ( outa == PDT::Colour0 && outb == PDT::Colour6bar ) {
return GeneralHardME::Colour3bar3barto16bar;
}
else if ( outa == PDT::Colour0 && outb == PDT::Colour3) {
return GeneralHardME::Colour3bar3barto13;
}
else if ( outb == PDT::Colour0 && outa == PDT::Colour3) {
return GeneralHardME::Colour3bar3barto31;
}
else if ( outa == PDT::Colour8 && outb == PDT::Colour3) {
return GeneralHardME::Colour3bar3barto83;
}
else if ( outb == PDT::Colour8 && outa == PDT::Colour3) {
return GeneralHardME::Colour3bar3barto38;
}
else
assert(false);
}
// incoming 3 3bar
else if(ina == PDT::Colour3 && inb == PDT::Colour3bar ) {
if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour33barto11;
}
else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
return GeneralHardME::Colour33barto33bar;
}
else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
return GeneralHardME::Colour33barto88;
}
else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour33barto81;
}
else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
return GeneralHardME::Colour33barto18;
}
else if( outa == PDT::Colour6 && outb == PDT::Colour6bar) {
return GeneralHardME::Colour33barto66bar;
}
else if( outa == PDT::Colour6bar && outb == PDT::Colour6) {
return GeneralHardME::Colour33barto6bar6;
}
else
assert(false);
}
// incoming 88
else if(ina == PDT::Colour8 && inb == PDT::Colour8 ) {
if( outa == PDT::Colour0 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour88to11;
}
else if( outa == PDT::Colour3 && outb == PDT::Colour3bar ) {
return GeneralHardME::Colour88to33bar;
}
else if( outa == PDT::Colour8 && outb == PDT::Colour8 ) {
return GeneralHardME::Colour88to88;
}
else if( outa == PDT::Colour8 && outb == PDT::Colour0 ) {
return GeneralHardME::Colour88to81;
}
else if( outa == PDT::Colour0 && outb == PDT::Colour8 ) {
return GeneralHardME::Colour88to18;
}
else if( outa == PDT::Colour6 && outb == PDT::Colour6bar ) {
return GeneralHardME::Colour88to66bar;
}
else
assert(false);
}
// incoming 38
else if(ina == PDT::Colour3 && inb == PDT::Colour8 ) {
if(outa == PDT::Colour3 && outb == PDT::Colour0) {
return GeneralHardME::Colour38to31;
}
else if(outa == PDT::Colour0 && outb == PDT::Colour3) {
return GeneralHardME::Colour38to13;
}
else if(outa == PDT::Colour3 && outb == PDT::Colour8) {
return GeneralHardME::Colour38to38;
}
else if(outa == PDT::Colour8 && outb == PDT::Colour3) {
return GeneralHardME::Colour38to83;
}
else if(outa == PDT::Colour3bar && outb == PDT::Colour6){
return GeneralHardME::Colour38to3bar6;
}
else if(outa == PDT::Colour6 && outb == PDT::Colour3bar) {
return GeneralHardME::Colour38to63bar;
}
else if(outa == PDT::Colour3bar && outb == PDT::Colour3bar) {
return GeneralHardME::Colour38to3bar3bar;
}
else
assert(false);
}
// incoming 3bar8
else if(ina == PDT::Colour3bar && inb == PDT::Colour8 ) {
if(outa == PDT::Colour3bar && outb == PDT::Colour0 ) {
return GeneralHardME::Colour3bar8to3bar1;
}
else if(outa == PDT::Colour0 && outb == PDT::Colour3bar) {
return GeneralHardME::Colour3bar8to13bar;
}
else if(outa == PDT::Colour3bar && outb == PDT::Colour8 ) {
return GeneralHardME::Colour3bar8to3bar8;
}
else if(outa == PDT::Colour8 && outb == PDT::Colour3bar) {
return GeneralHardME::Colour3bar8to83bar;
}
else if(outa == PDT::Colour3 && outb == PDT::Colour3) {
return GeneralHardME::Colour3bar8to33;
}
else
assert(false);
}
// unknown colour flow
else
assert(false);
return GeneralHardME::UNDEFINED;
}
void HardProcessConstructor::fixFSOrder(HPDiagram & diag) {
tcPDPtr psa = getParticleData(diag.incoming.first);
tcPDPtr psb = getParticleData(diag.incoming.second);
tcPDPtr psc = getParticleData(diag.outgoing.first);
tcPDPtr psd = getParticleData(diag.outgoing.second);
//fix a spin order
if( psc->iSpin() < psd->iSpin() ) {
swap(diag.outgoing.first, diag.outgoing.second);
if(diag.channelType == HPDiagram::tChannel) {
diag.ordered.second = !diag.ordered.second;
}
return;
}
if( psc->iSpin() == psd->iSpin() &&
psc->id() < 0 && psd->id() > 0 ) {
swap(diag.outgoing.first, diag.outgoing.second);
if(diag.channelType == HPDiagram::tChannel) {
diag.ordered.second = !diag.ordered.second;
}
return;
}
}
void HardProcessConstructor::assignToCF(HPDiagram & diag) {
if(diag.channelType == HPDiagram::tChannel) {
if(diag.ordered.second) tChannelCF(diag);
else uChannelCF(diag);
}
else if(diag.channelType == HPDiagram::sChannel) {
sChannelCF(diag);
}
else if (diag.channelType == HPDiagram::fourPoint) {
fourPointCF(diag);
}
else
assert(false);
}
void HardProcessConstructor::tChannelCF(HPDiagram & diag) {
tcPDPtr ia = getParticleData(diag.incoming.first );
tcPDPtr ib = getParticleData(diag.incoming.second);
tcPDPtr oa = getParticleData(diag.outgoing.first );
tcPDPtr ob = getParticleData(diag.outgoing.second);
PDT::Colour ina = ia->iColour();
PDT::Colour inb = ib->iColour();
PDT::Colour outa = oa->iColour();
PDT::Colour outb = ob->iColour();
vector<CFPair> cfv(1, make_pair(0, 1.));
if(diag.intermediate->iColour() == PDT::Colour0) {
if(ina==PDT::Colour0) {
cfv[0] = make_pair(0, 1);
}
else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
if( inb == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
cfv[0] = make_pair(2, 1);
}
else if(inb==PDT::Colour8) {
cfv[0] = make_pair(2, 1);
}
}
else if(ina==PDT::Colour8) {
if( inb == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
cfv[0] = make_pair(2, 1);
}
else if(inb==PDT::Colour8) {
cfv[0] = make_pair(7, 1);
}
}
}
else if(diag.intermediate->iColour() == PDT::Colour8) {
if(ina==PDT::Colour8&&outa==PDT::Colour8&&
inb==PDT::Colour8&&outb==PDT::Colour8) {
cfv[0]=make_pair(2, 2.);
cfv.push_back(make_pair(3, -2.));
cfv.push_back(make_pair(1, -2.));
cfv.push_back(make_pair(4, 2.));
}
else if(ina==PDT::Colour8&&outa==PDT::Colour0&&
inb==PDT::Colour8&&outb==PDT::Colour8&&
oa->iSpin()==PDT::Spin0) {
cfv[0] = make_pair(0,-1);
}
else if(ina==PDT::Colour8&&outa==PDT::Colour8&&
inb==PDT::Colour8&&outb==PDT::Colour0&&
ob->iSpin()==PDT::Spin0) {
cfv[0] = make_pair(0,-1);
}
}
else if(diag.intermediate->iColour() == PDT::Colour3 ||
diag.intermediate->iColour() == PDT::Colour3bar) {
if(outa == PDT::Colour0 || outb == PDT::Colour0) {
if( outa != PDT::Colour6 && outb != PDT::Colour6 &&
outa != PDT::Colour6bar && outb != PDT::Colour6bar) {
cfv[0] = make_pair(0,1.);
}
else {
cfv[0] = make_pair(0,0.5);
cfv.push_back(make_pair(1,0.5));
}
}
else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
cfv[0] = make_pair(4,1.);
cfv.push_back(make_pair(5,1.));
}
else if(outa==PDT::Colour6 && outb==PDT::Colour6bar) {
cfv[0] = make_pair(4, 1.);
for(unsigned int ix=5;ix<8;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(outa==PDT::Colour6 || outa ==PDT::Colour6bar ||
outb==PDT::Colour6 || outb ==PDT::Colour6bar ) {
assert(false);
}
else if(ina==PDT::Colour3 && inb==PDT::Colour3 ) {
if((outa==PDT::Colour0 && outb==PDT::Colour3bar)||
(outb==PDT::Colour0 && outa==PDT::Colour3bar))
cfv[0] = make_pair(0,1.);
else if((outa==PDT::Colour8 && outb==PDT::Colour3bar)||
(outb==PDT::Colour8 && outa==PDT::Colour3bar))
cfv[0] = make_pair(1,1.);
}
else if(ina==PDT::Colour3bar && inb==PDT::Colour3bar ) {
if((outa==PDT::Colour0 && outb==PDT::Colour3)||
(outb==PDT::Colour0 && outa==PDT::Colour3))
cfv[0] = make_pair(0,1.);
else if((outa==PDT::Colour8 && outb==PDT::Colour3)||
(outb==PDT::Colour8 && outa==PDT::Colour3))
cfv[0] = make_pair(1,1.);
}
else if((ina==PDT::Colour3 && inb==PDT::Colour8) ||
(ina==PDT::Colour3bar && inb==PDT::Colour8) ||
(inb==PDT::Colour3 && ina==PDT::Colour8) ||
(inb==PDT::Colour3bar && ina==PDT::Colour8) ) {
if((outa==PDT::Colour3 && outb==PDT::Colour3 ) ||
(outa==PDT::Colour3bar && outb==PDT::Colour3bar)) {
cfv[0] = make_pair(1,1.);
}
}
}
else if(diag.intermediate->iColour() == PDT::Colour6 ||
diag.intermediate->iColour() == PDT::Colour6bar) {
if(ina==PDT::Colour8 && inb==PDT::Colour8) {
cfv[0] = make_pair(0, 1.);
for(unsigned int ix=1;ix<4;++ix)
cfv.push_back(make_pair(ix,1.));
for(unsigned int ix=4;ix<8;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
cfv[0] = make_pair(0,1.);
for(unsigned int ix=1;ix<4;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
cfv[0] = make_pair(4,1.);
cfv.push_back(make_pair(5,1.));
}
}
diag.colourFlow = cfv;
}
void HardProcessConstructor::uChannelCF(HPDiagram & diag) {
PDT::Colour offshell = diag.intermediate->iColour();
PDT::Colour ina = getParticleData(diag.incoming.first )->iColour();
PDT::Colour inb = getParticleData(diag.incoming.second)->iColour();
PDT::Colour outa = getParticleData(diag.outgoing.first )->iColour();
PDT::Colour outb = getParticleData(diag.outgoing.second)->iColour();
vector<CFPair> cfv(1, make_pair(1, 1.));
if(offshell == PDT::Colour8) {
if(outa == PDT::Colour0 &&
outb == PDT::Colour0) {
cfv[0].first = 0;
}
else if( outa != outb ) {
if(outa == PDT::Colour0 ||
outb == PDT::Colour0) {
cfv[0].first = 0;
}
else if(ina == PDT::Colour3 && inb == PDT::Colour8 &&
outb == PDT::Colour3 && outa == PDT::Colour8) {
tPDPtr off = diag.intermediate;
if(off->CC()) off=off->CC();
if(off->iSpin()!=PDT::Spin1Half ||
diag.vertices.second->allowed(off->id(),diag.outgoing.first,diag.incoming.second)) {
cfv[0].first = 0;
cfv.push_back(make_pair(1, -1.));
}
else {
cfv[0].first = 1;
cfv.push_back(make_pair(0, -1.));
}
}
else if(ina == PDT::Colour3bar && inb == PDT::Colour8 &&
outb == PDT::Colour3bar && outa == PDT::Colour8) {
tPDPtr off = diag.intermediate;
if(off->CC()) off=off->CC();
if(off->iSpin()!=PDT::Spin1Half ||
diag.vertices.second->allowed(diag.outgoing.first,off->id(),diag.incoming.second)) {
cfv[0].first = 0;
cfv.push_back(make_pair(1, -1.));
}
else {
cfv[0].first = 1;
cfv.push_back(make_pair(0, -1.));
}
}
else {
cfv[0].first = 0;
cfv.push_back(make_pair(1, -1.));
}
}
else if(outa==PDT::Colour8&&ina==PDT::Colour8) {
cfv[0]=make_pair(4, 2.);
cfv.push_back(make_pair(5, -2.));
cfv.push_back(make_pair(0, -2.));
cfv.push_back(make_pair(2, 2.));
}
}
else if(offshell == PDT::Colour3 || offshell == PDT::Colour3bar) {
if( outa == PDT::Colour0 || outb == PDT::Colour0 ) {
if( outa != PDT::Colour6 && outb != PDT::Colour6 &&
outa != PDT::Colour6bar && outb != PDT::Colour6bar) {
cfv[0] = make_pair(0,1.);
}
else {
cfv[0] = make_pair(0,0.5);
cfv.push_back(make_pair(1,0.5));
}
}
else if(outa==PDT::Colour3bar && outb==PDT::Colour6) {
cfv[0] = make_pair(4,1.);
cfv.push_back(make_pair(5,1.));
}
else if(outa==PDT::Colour6 && outb==PDT::Colour3bar) {
cfv[0] = make_pair(0,1.);
for(int ix=0; ix<4;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(outa==PDT::Colour6bar && outb==PDT::Colour6) {
cfv[0] = make_pair(4,1.);
for(int ix=5; ix<8;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(ina==PDT::Colour0 && inb==PDT::Colour0) {
cfv[0] = make_pair(0,1.);
}
else if(ina==PDT::Colour3 && inb==PDT::Colour3 ) {
if((outa==PDT::Colour0 && outb==PDT::Colour3bar)||
(outb==PDT::Colour0 && outa==PDT::Colour3bar))
cfv[0] = make_pair(0,1.);
else if((outa==PDT::Colour8 && outb==PDT::Colour3bar)||
(outb==PDT::Colour8 && outa==PDT::Colour3bar))
cfv[0] = make_pair(2,1.);
}
else if(ina==PDT::Colour3bar && inb==PDT::Colour3bar ) {
if((outa==PDT::Colour0 && outb==PDT::Colour3)||
(outb==PDT::Colour0 && outa==PDT::Colour3))
cfv[0] = make_pair(0,1.);
else if((outa==PDT::Colour8 && outb==PDT::Colour3)||
(outb==PDT::Colour8 && outa==PDT::Colour3))
cfv[0] = make_pair(2,1.);
}
else if(((ina==PDT::Colour3 && inb==PDT::Colour8) ||
(ina==PDT::Colour3bar && inb==PDT::Colour8) ||
(inb==PDT::Colour3 && ina==PDT::Colour8) ||
(inb==PDT::Colour3bar && ina==PDT::Colour8)) &&
((outa==PDT::Colour3 && outb==PDT::Colour3 ) ||
(outa==PDT::Colour3bar && outb==PDT::Colour3bar))) {
cfv[0] = make_pair(2, 1.);
}
else if(( ina==PDT::Colour3 && inb==PDT::Colour3bar &&
outa==PDT::Colour3 && outb==PDT::Colour3bar)) {
cfv[0] = make_pair(2, 1.);
cfv.push_back(make_pair(3,-1.));
}
}
else if( offshell == PDT::Colour0 ) {
if(ina==PDT::Colour0) {
cfv[0] = make_pair(0, 1);
}
else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
if( inb == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
cfv[0] = make_pair(3, 1);
}
else if(inb==PDT::Colour8) {
cfv[0] = make_pair(2, 1);
}
}
else if(ina==PDT::Colour8) {
if( inb == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(inb==PDT::Colour3 || outb==PDT::Colour3bar) {
cfv[0] = make_pair(2, 1);
}
else if(inb==PDT::Colour8) {
cfv[0] = make_pair(8, 1);
}
}
}
else if(diag.intermediate->iColour() == PDT::Colour6 ||
diag.intermediate->iColour() == PDT::Colour6bar) {
if(ina==PDT::Colour8 && inb==PDT::Colour8) {
cfv[0] = make_pair(0, 1.);
for(unsigned int ix=1;ix<4;++ix)
cfv.push_back(make_pair(ix,1.));
for(unsigned int ix=8;ix<12;++ix)
cfv.push_back(make_pair(ix,1.));
}
else if(outa==PDT::Colour3bar && outa==PDT::Colour6) {
cfv[0] = make_pair(4, 1.);
cfv.push_back(make_pair(5,1.));
}
else if(outa==PDT::Colour6 && outa==PDT::Colour3bar) {
cfv[0] = make_pair(0, 1.);
for(unsigned int ix=1;ix<4;++ix)
cfv.push_back(make_pair(ix,1.));
}
}
diag.colourFlow = cfv;
}
void HardProcessConstructor::sChannelCF(HPDiagram & diag) {
tcPDPtr pa = getParticleData(diag.incoming.first);
tcPDPtr pb = getParticleData(diag.incoming.second);
PDT::Colour ina = pa->iColour();
PDT::Colour inb = pb->iColour();
PDT::Colour offshell = diag.intermediate->iColour();
tcPDPtr pc = getParticleData(diag.outgoing.first);
tcPDPtr pd = getParticleData(diag.outgoing.second);
PDT::Colour outa = pc->iColour();
PDT::Colour outb = pd->iColour();
vector<CFPair> cfv(1);
if(offshell == PDT::Colour8) {
if(ina == PDT::Colour0 || inb == PDT::Colour0 ||
outa == PDT::Colour0 || outb == PDT::Colour0) {
cfv[0] = make_pair(0, 1);
}
else {
bool incol = ina == PDT::Colour8 && inb == PDT::Colour8;
bool outcol = outa == PDT::Colour8 && outb == PDT::Colour8;
bool intrip = ina == PDT::Colour3 && inb == PDT::Colour3bar;
bool outtrip = outa == PDT::Colour3 && outb == PDT::Colour3bar;
bool outsex = outa == PDT::Colour6 && outb == PDT::Colour6bar;
bool outsexb = outa == PDT::Colour6bar && outb == PDT::Colour6;
if(incol || outcol) {
// Require an additional minus sign for a scalar/fermion
// 33bar final state due to the way the vertex rules are defined.
int prefact(1);
if( ((pc->iSpin() == PDT::Spin1Half && pd->iSpin() == PDT::Spin1Half) ||
(pc->iSpin() == PDT::Spin0 && pd->iSpin() == PDT::Spin0 )) &&
(outa == PDT::Colour3 && outb == PDT::Colour3bar) )
prefact = -1;
if(incol && outcol) {
cfv[0] = make_pair(0, -2.);
cfv.push_back(make_pair(1, 2.));
cfv.push_back(make_pair(3, 2.));
cfv.push_back(make_pair(5, -2.));
}
else if(incol && outsex) {
cfv[0].first = 4;
cfv[0].second = prefact;
for(unsigned int ix=1;ix<4;++ix)
cfv.push_back(make_pair(4+ix, prefact));
for(unsigned int ix=0;ix<4;++ix)
cfv.push_back(make_pair(8+ix,-prefact));
}
else {
cfv[0].first = 0;
cfv[0].second = -prefact;
cfv.push_back(make_pair(1, prefact));
}
}
else if( ( intrip && !outtrip ) ||
( !intrip && outtrip ) ) {
if(!outsex)
cfv[0] = make_pair(0, 1);
else {
cfv[0] = make_pair(0, 1.);
for(unsigned int ix=0;ix<3;++ix)
cfv.push_back(make_pair(ix+1, 1.));
}
}
else if((intrip && outsex) || (intrip && outsexb)) {
cfv[0] = make_pair(0,1.);
for(int ix=1; ix<4; ++ix)
cfv.push_back(make_pair(ix,1.));
}
else
cfv[0] = make_pair(1, 1);
}
}
else if(offshell == PDT::Colour0) {
if( ina == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(ina==PDT::Colour3 || ina==PDT::Colour3bar) {
if( outa == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(outa==PDT::Colour3 || outa==PDT::Colour3bar) {
cfv[0] = make_pair(3, 1);
}
else if(outa==PDT::Colour8) {
cfv[0] = make_pair(2, 1);
}
else if(outa==PDT::Colour6 || outa==PDT::Colour6bar) {
cfv[0] = make_pair(8, 1.);
cfv.push_back(make_pair(9,1.));
}
else
assert(false);
}
else if(ina==PDT::Colour8) {
if( outa == PDT::Colour0 ) {
cfv[0] = make_pair(0, 1);
}
else if(outa==PDT::Colour3 || outb==PDT::Colour3bar) {
cfv[0] = make_pair(2, 1);
}
else if(outa==PDT::Colour8) {
cfv[0] = make_pair(6, 1);
}
}
}
else if(offshell == PDT::Colour3 || offshell == PDT::Colour3bar) {
if(outa == PDT::Colour6 || outa == PDT::Colour6bar ||
outb == PDT::Colour6bar || outb == PDT::Colour6) {
cfv[0] = make_pair(6, 1.);
cfv.push_back(make_pair(7,1.));
}
else if((ina == PDT::Colour3 && inb == PDT::Colour3) ||
(ina == PDT::Colour3bar && inb == PDT::Colour3bar)) {
if((outa == PDT::Colour3 && outb == PDT::Colour3 ) ||
(outa == PDT::Colour3bar && outb == PDT::Colour3bar)) {
cfv[0] = make_pair(2, 1.);
cfv.push_back(make_pair(3,-1.));
}
else
cfv[0] = make_pair(0,1.);
}
else if(((ina==PDT::Colour3 && inb==PDT::Colour8) ||
(ina==PDT::Colour3bar && inb==PDT::Colour8) ||
(inb==PDT::Colour3 && ina==PDT::Colour8) ||
(inb==PDT::Colour3bar && ina==PDT::Colour8) ) &&
((outa==PDT::Colour3 && outb==PDT::Colour3 ) ||
(outa==PDT::Colour3bar && outb==PDT::Colour3bar))) {
cfv[0] = make_pair(0,1.);
}
else {
if(outa == PDT::Colour0 || outb == PDT::Colour0)
cfv[0] = make_pair(0, 1);
else
cfv[0] = make_pair(1, 1);
}
}
else if( offshell == PDT::Colour6 || offshell == PDT::Colour6bar) {
if((ina == PDT::Colour3 && inb == PDT::Colour3 &&
outa == PDT::Colour3 && outb == PDT::Colour3 ) ||
(ina == PDT::Colour3bar && inb == PDT::Colour3bar &&
outa == PDT::Colour3bar && outb == PDT::Colour3bar)) {
cfv[0] = make_pair(2,0.5);
cfv.push_back(make_pair(3,0.5));
}
else if((ina == PDT::Colour3 && inb == PDT::Colour3 &&
((outa == PDT::Colour6 && outb == PDT::Colour0)||
(outb == PDT::Colour6 && outa == PDT::Colour0))) ||
(ina == PDT::Colour3bar && inb == PDT::Colour3bar &&
((outa == PDT::Colour6bar && outb == PDT::Colour0)||
(outb == PDT::Colour6bar && outa == PDT::Colour0)))) {
cfv[0] = make_pair(0,0.5);
cfv.push_back(make_pair(1,0.5));
}
else
assert(false);
}
else {
if(outa == PDT::Colour0 || outb == PDT::Colour0)
cfv[0] = make_pair(0, 1);
else
cfv[0] = make_pair(1, 1);
}
diag.colourFlow = cfv;
}
void HardProcessConstructor::fourPointCF(HPDiagram & diag) {
// count the colours
unsigned int noct(0),ntri(0),nsng(0),nsex(0),nf(0);
for(unsigned int ix=0;ix<4;++ix) {
tcPDPtr pd = getParticleData(diag.ids[ix]);
PDT::Colour col = pd->iColour();
if(col==PDT::Colour0) ++nsng;
else if(col==PDT::Colour3||col==PDT::Colour3bar) ++ntri;
else if(col==PDT::Colour8) ++noct;
else if(col==PDT::Colour6||col==PDT::Colour6bar) ++nsex;
if(pd->iSpin()==2) nf+=1;
}
if(nsng==4 || (ntri==2&&nsng==2) ||
(noct==3 && nsng==1) ||
(ntri==2 && noct==1 && nsng==1) ||
(noct == 2 && nsng == 2) ) {
vector<CFPair> cfv(1,make_pair(0,1));
diag.colourFlow = cfv;
}
else if(noct==4) {
// flows for SSVV, VVVV is handled in me class
vector<CFPair> cfv(6);
cfv[0] = make_pair(0, -2.);
cfv[1] = make_pair(1, -2.);
cfv[2] = make_pair(2, +4.);
cfv[3] = make_pair(3, -2.);
cfv[4] = make_pair(4, +4.);
cfv[5] = make_pair(5, -2.);
diag.colourFlow = cfv;
}
else if(ntri==2&&noct==2) {
vector<CFPair> cfv(2);
cfv[0] = make_pair(0, 1);
cfv[1] = make_pair(1, 1);
if(nf==2) cfv[1].second = -1.;
diag.colourFlow = cfv;
}
else if(nsex==2&&noct==2) {
vector<CFPair> cfv;
for(unsigned int ix=0;ix<4;++ix)
cfv.push_back(make_pair(ix ,2.));
for(unsigned int ix=0;ix<8;++ix)
cfv.push_back(make_pair(4+ix,1.));
diag.colourFlow = cfv;
}
else
assert(false);
}
namespace {
// Helper functor for find_if in duplicate function.
class SameDiagramAs {
public:
SameDiagramAs(const HPDiagram & diag) : a(diag) {}
bool operator()(const HPDiagram & b) const {
return a == b;
}
private:
HPDiagram a;
};
}
bool HardProcessConstructor::duplicate(const HPDiagram & diag,
const HPDVector & group) const {
//find if a duplicate diagram exists
HPDVector::const_iterator it =
find_if(group.begin(), group.end(), SameDiagramAs(diag));
return it != group.end();
}
bool HardProcessConstructor::checkOrder(const HPDiagram & diag) const {
for(map<string,pair<unsigned int,unsigned int> >::const_iterator it=model_->couplings().begin();
it!=model_->couplings().end();++it) {
- cerr << "testing in check order loop? " << __LINE__ << "\n";
unsigned int order=0;
- cerr << "testing in check order loop? " << __LINE__ << " " << it->second.first << " "
- << it->first << "\n";
if(diag.vertices.first ) order += diag.vertices.first ->orderInCoupling(it->second.first);
- cerr << "testing in check order loop? " << __LINE__ << "\n";
if(diag.vertices.second&&diag.vertices.first->getNpoint()==3)
order += diag.vertices.second->orderInCoupling(it->second.first);
- cerr << "testing in check order loop? " << __LINE__ << "\n";
if(order>it->second.second) return false;
- cerr << "testing in check order loop? " << __LINE__ << "\n";
}
return true;
}
diff --git a/Models/General/ModelGenerator.cc b/Models/General/ModelGenerator.cc
--- a/Models/General/ModelGenerator.cc
+++ b/Models/General/ModelGenerator.cc
@@ -1,566 +1,555 @@
// -*- C++ -*-
//
// ModelGenerator.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the ModelGenerator class.
//
#include "ModelGenerator.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_;
}
// 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();
- cerr << "testing in doinit " << __LINE__ << "\n";
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_);
}
- cerr << "testing in doinit " << __LINE__ << "\n";
// 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());
}
- cerr << "testing in doinit " << __LINE__ << "\n";
-
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();
// Now switch off the modes if needed
for(DecaySet::const_iterator it=parent->decayModes().begin();
it!=parent->decayModes().end();++it) {
if( _theDecayConstructor->disableDecayMode((**it).tag()) ) {
DMPtr mode = *it;
generator()->preinitInterface(mode, "Active", "set", "No");
if(mode->CC()) {
DMPtr CCmode = mode->CC();
generator()->preinitInterface(CCmode, "Active", "set", "No");
}
}
}
- cerr << "testing in doinit " << __LINE__ << "\n";
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());
}
}
}
- cerr << "testing in doinit " << __LINE__ << "\n";
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();
}
- cerr << "testing in doinit " << __LINE__ << "\n";
// loop again to initialise mass and width generators
// switch off modes and write output
for(PDVector::iterator pit = particles_.begin();
pit != particles_.end(); ++pit) {
tPDPtr parent = *pit;
if(parent->widthGenerator())
parent->widthGenerator()->init();
if(parent->massGenerator())
parent->massGenerator()->init();
// output the modes if needed
if( !parent->decaySelector().empty() ) {
if ( decayOutput_ == 2 )
writeDecayModes(ofs, parent);
else
writeDecayModes(os, parent);
}
}
- cerr << "testing in doinit " << __LINE__ << "\n";
//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();
- cerr << "before construct " << ix << " " << hardProcessConstructors_[ix]->fullName() << "\n";
hardProcessConstructors_[ix]->constructDiagrams();
- cerr << "before construct " << ix << "\n";
}
- cerr << "testing in doinit " << __LINE__ << "\n";
}
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, "Active", "set", "No");
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 no modes left set stable
if(newwidth==ZERO) {
parent->stable(true);
parent->width(ZERO);
parent->widthCut(ZERO);
parent->massGenerator(tGenericMassGeneratorPtr());
parent->widthGenerator(tGenericWidthGeneratorPtr());
}
// otherwise rescale if needed
else 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(parent->decayModes().begin(),
parent->decayModes().end());
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"
<< std::left << std::setw(20) << "BR" << "Yes/No\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
<< std::left << std::setw(20) << (**dit).brat()
<< ((**dit).on() ? "Yes" : "No" ) << '\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/TwoToTwoProcessConstructor.cc b/Models/General/TwoToTwoProcessConstructor.cc
--- a/Models/General/TwoToTwoProcessConstructor.cc
+++ b/Models/General/TwoToTwoProcessConstructor.cc
@@ -1,666 +1,643 @@
// -*- C++ -*-
//
// TwoToTwoProcessConstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2017 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
//
// This is the implementation of the non-inlined, non-templated member
// functions of the TwoToTwoProcessConstructor class.
//
#include "TwoToTwoProcessConstructor.h"
#include "ThePEG/Utilities/DescribeClass.h"
#include "ThePEG/Interface/ClassDocumentation.h"
#include "ThePEG/Persistency/PersistentOStream.h"
#include "ThePEG/Persistency/PersistentIStream.h"
#include "ThePEG/Interface/RefVector.h"
#include "ThePEG/Interface/Reference.h"
#include "ThePEG/Interface/Parameter.h"
#include "ThePEG/Interface/Switch.h"
#include <sstream>
using std::stringstream;
using namespace Herwig;
TwoToTwoProcessConstructor::TwoToTwoProcessConstructor() :
Nout_(0), nv_(0), allDiagrams_(true),
processOption_(0), scaleChoice_(0), scaleFactor_(1.)
{}
IBPtr TwoToTwoProcessConstructor::clone() const {
return new_ptr(*this);
}
IBPtr TwoToTwoProcessConstructor::fullclone() const {
return new_ptr(*this);
}
void TwoToTwoProcessConstructor::doinit() {
HardProcessConstructor::doinit();
if((processOption_==2 || processOption_==4) &&
outgoing_.size()!=2)
throw InitException()
<< "Exclusive processes require exactly"
<< " two outgoing particles but " << outgoing_.size()
<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
<< Exception::runerror;
if(processOption_==4 && incoming_.size()!=2)
throw InitException()
<< "Exclusive processes require exactly"
<< " two incoming particles but " << incoming_.size()
<< "have been inserted in TwoToTwoProcessConstructor::doinit()."
<< Exception::runerror;
Nout_ = outgoing_.size();
PDVector::size_type ninc = incoming_.size();
// exit if nothing to do
if(Nout_==0||ninc==0) return;
//create vector of initial-state pairs
for(PDVector::size_type i = 0; i < ninc; ++i) {
for(PDVector::size_type j = 0; j < ninc; ++j) {
tPDPair inc = make_pair(incoming_[i], incoming_[j]);
if( (inc.first->iSpin() > inc.second->iSpin()) ||
(inc.first->iSpin() == inc.second->iSpin() &&
inc.first->id() < inc.second->id()) )
swap(inc.first, inc.second);
if( !HPC_helper::duplicateIncoming(inc,incPairs_) ) {
incPairs_.push_back(inc);
}
}
}
// excluded vertices
excludedVertexSet_ =
set<VertexBasePtr>(excludedVertexVector_.begin(),
excludedVertexVector_.end());
}
void TwoToTwoProcessConstructor::persistentOutput(PersistentOStream & os) const {
os << vertices_ << incoming_ << outgoing_
<< allDiagrams_ << processOption_
<< scaleChoice_ << scaleFactor_ << excluded_ << excludedExternal_
<< excludedVertexVector_ << excludedVertexSet_;
}
void TwoToTwoProcessConstructor::persistentInput(PersistentIStream & is, int) {
is >> vertices_ >> incoming_ >> outgoing_
>> allDiagrams_ >> processOption_
>> scaleChoice_ >> scaleFactor_ >> excluded_ >> excludedExternal_
>> excludedVertexVector_ >> excludedVertexSet_;
}
// The following static variable is needed for the type
// description system in ThePEG.
DescribeClass<TwoToTwoProcessConstructor,HardProcessConstructor>
describeHerwigTwoToTwoProcessConstructor("Herwig::TwoToTwoProcessConstructor", "Herwig.so");
void TwoToTwoProcessConstructor::Init() {
static ClassDocumentation<TwoToTwoProcessConstructor> documentation
("TwoToTwoProcessConstructor constructs the possible diagrams for "
"a process given the external particles");
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceIn
("Incoming",
"Pointers to incoming particles",
&TwoToTwoProcessConstructor::incoming_, -1, false, false, true, false);
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceOut
("Outgoing",
"Pointers to incoming particles",
&TwoToTwoProcessConstructor::outgoing_, -1, false, false, true, false);
static Switch<TwoToTwoProcessConstructor,bool> interfaceIncludeAllDiagrams
("IncludeEW",
"Switch to decide which diagrams to include in ME calc.",
&TwoToTwoProcessConstructor::allDiagrams_, true, false, false);
static SwitchOption interfaceIncludeAllDiagramsNo
(interfaceIncludeAllDiagrams,
"No",
"Only include QCD diagrams",
false);
static SwitchOption interfaceIncludeAllDiagramsYes
(interfaceIncludeAllDiagrams,
"Yes",
"Include EW+QCD.",
true);
static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceProcesses
("Processes",
"Whether to generate inclusive or exclusive processes",
&TwoToTwoProcessConstructor::processOption_, 0, false, false);
static SwitchOption interfaceProcessesSingleParticleInclusive
(interfaceProcesses,
"SingleParticleInclusive",
"Require at least one particle from the list of outgoing particles"
" in the hard process",
0);
static SwitchOption interfaceProcessesTwoParticleInclusive
(interfaceProcesses,
"TwoParticleInclusive",
"Require that both the particles in the hard processes are in the"
" list of outgoing particles",
1);
static SwitchOption interfaceProcessesExclusive
(interfaceProcesses,
"Exclusive",
"Require that both the particles in the hard processes are in the"
" list of outgoing particles in every hard process",
2);
static SwitchOption interfaceProcessesVeryExclusive
(interfaceProcesses,
"VeryExclusive",
"Require that both the incoming and outgoing particles in the hard processes are in the"
" list of outgoing particles in every hard process",
4);
static Switch<TwoToTwoProcessConstructor,unsigned int> interfaceScaleChoice
("ScaleChoice",
"&TwoToTwoProcessConstructor::scaleChoice_",
&TwoToTwoProcessConstructor::scaleChoice_, 0, false, false);
static SwitchOption interfaceScaleChoiceDefault
(interfaceScaleChoice,
"Default",
"Use if sHat if intermediates all colour neutral, otherwise the transverse mass",
0);
static SwitchOption interfaceScaleChoicesHat
(interfaceScaleChoice,
"sHat",
"Always use sHat",
1);
static SwitchOption interfaceScaleChoiceTransverseMass
(interfaceScaleChoice,
"TransverseMass",
"Always use the transverse mass",
2);
static SwitchOption interfaceScaleChoiceGeometicMean
(interfaceScaleChoice,
"MaxMT",
"Use the maximum of m^2+p_T^2 for the two particles",
3);
static Parameter<TwoToTwoProcessConstructor,double> interfaceScaleFactor
("ScaleFactor",
"The prefactor used in the scale calculation. The scale used is"
" that defined by scaleChoice multiplied by this prefactor",
&TwoToTwoProcessConstructor::scaleFactor_, 1.0, 0.0, 10.0,
false, false, Interface::limited);
static RefVector<TwoToTwoProcessConstructor,ThePEG::ParticleData> interfaceExcluded
("Excluded",
"Particles which are not allowed as intermediates",
&TwoToTwoProcessConstructor::excluded_, -1, false, false, true, false, false);
static RefVector<TwoToTwoProcessConstructor,ParticleData> interfaceExcludedExternal
("ExcludedExternal",
"Particles which are not allowed as outgoing particles",
&TwoToTwoProcessConstructor::excludedExternal_, -1,
false, false, true, false, false);
static RefVector<TwoToTwoProcessConstructor,VertexBase> interfaceExcludedVertices
("ExcludedVertices",
"Vertices which are not included in the 2 -> 2 scatterings",
&TwoToTwoProcessConstructor::excludedVertexVector_, -1, false, false, true, true, false);
}
void TwoToTwoProcessConstructor::constructDiagrams() {
- cerr << "testing in construct " << __LINE__ << "\n";
if(incPairs_.empty() || outgoing_.empty() || !subProcess() ) return;
nv_ = model()->numberOfVertices();
//make sure vertices are initialised
for(unsigned int ix = 0; ix < nv_; ++ix ) {
VertexBasePtr vertex = model()->vertex(ix);
if(excludedVertexSet_.find(vertex) !=
excludedVertexSet_.end()) continue;
vertices_.push_back(vertex);
}
nv_ = vertices_.size();
- cerr << "testing in construct " << __LINE__ << "\n";
//Create necessary diagrams
vector<tcPDPair>::size_type is;
PDVector::size_type os;
- cerr << "testing in construct " << __LINE__ << "\n";
for(is = 0; is < incPairs_.size(); ++is) {
tPDPair ppi = incPairs_[is];
for(os = 0; os < Nout_; ++os) {
long fs = outgoing_[os]->id();
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexA = vertices_[iv];
//This skips an effective vertex and the EW ones if
// we only want the strong diagrams
if( !allDiagrams_ && vertexA->orderInGs() == 0 )
continue;
if(vertexA->getNpoint() == 3) {
//scattering diagrams
- cerr << "testing before t\n";
createTChannels(ppi, fs, vertexA);
- cerr << "testing after t\n";
//resonance diagrams
if( vertexA->isIncoming(ppi.first) &&
vertexA->isIncoming(ppi.second) )
createSChannels(ppi, fs, vertexA);
- cerr << "testing after s\n";
}
else {
- cerr << "testing before 4\n";
makeFourPointDiagrams(ppi.first->id(), ppi.second->id(),
fs, vertexA);
- cerr << "testing after 4\n";
}
}
}
}
- cerr << "testing in construct " << __LINE__ << "\n";
// need to find all of the diagrams that relate to the same process
// first insert them into a map which uses the '<' operator
// to sort the diagrams
multiset<HPDiagram> grouped;
HPDVector::iterator dit = processes_.begin();
HPDVector::iterator dend = processes_.end();
for( ; dit != dend; ++dit) {
grouped.insert(*dit);
}
- cerr << "testing in construct " << __LINE__ << "\n";
assert( processes_.size() == grouped.size() );
processes_.clear();
typedef multiset<HPDiagram>::const_iterator set_iter;
set_iter it = grouped.begin(), iend = grouped.end();
- cerr << "testing in construct " << __LINE__ << "\n";
while( it != iend ) {
pair<set_iter,set_iter> range = grouped.equal_range(*it);
set_iter itb = range.first;
HPDVector process;
for( ; itb != range.second; ++itb ) {
process.push_back(*itb);
}
// if inclusive enforce the exclusivity
if(processOption_==2 || processOption_==4) {
if(!((process[0].outgoing. first==outgoing_[0]->id()&&
process[0].outgoing.second==outgoing_[1]->id())||
(process[0].outgoing. first==outgoing_[1]->id()&&
process[0].outgoing.second==outgoing_[0]->id()))) {
process.clear();
it = range.second;
continue;
}
if(processOption_==4) {
if(!((process[0].incoming. first==incoming_[0]->id()&&
process[0].incoming.second==incoming_[1]->id())||
(process[0].incoming. first==incoming_[1]->id()&&
process[0].incoming.second==incoming_[0]->id()))) {
process.clear();
it = range.second;
continue;
}
}
}
// check no zero width s-channel intermediates
for( dit=process.begin(); dit != process.end(); ++dit) {
tPDPtr out1 = getParticleData(dit->outgoing.first );
tPDPtr out2 = getParticleData(dit->outgoing.second);
if(dit->channelType == HPDiagram::sChannel &&
dit->intermediate->width()==ZERO &&
dit->intermediate->mass() > out1->mass()+ out2->mass()) {
tPDPtr in1 = getParticleData(dit->incoming.first );
tPDPtr in2 = getParticleData(dit->incoming.second);
throw Exception() << "Process with zero width resonant intermediates\n"
<< dit->intermediate->PDGName()
<< " can be on-shell in the process "
<< in1 ->PDGName() << " " << in2->PDGName() << " -> "
<< out1->PDGName() << " " << out2->PDGName()
<< " but has zero width.\nEither set the width, enable "
<< "calculation of its decays, and hence the width,\n"
<< "or disable it as a potential intermediate using\n"
<< "insert " << fullName() << ":Excluded 0 "
<< dit->intermediate->fullName() << "\n---\n"
<< Exception::runerror;
}
}
if(find(excludedExternal_.begin(),excludedExternal_.end(),
getParticleData(process[0].outgoing. first))!=excludedExternal_.end()) {
process.clear();
it = range.second;
continue;
}
if(find(excludedExternal_.begin(),excludedExternal_.end(),
getParticleData(process[0].outgoing.second))!=excludedExternal_.end()) {
process.clear();
it = range.second;
continue;
}
// finally if the process is allow assign the colour flows
for(unsigned int ix=0;ix<process.size();++ix) assignToCF(process[ix]);
// create the matrix element
createMatrixElement(process);
process.clear();
it = range.second;
}
- cerr << "testing in construct " << __LINE__ << "\n";
}
void TwoToTwoProcessConstructor::
createSChannels(tcPDPair inpp, long fs, tVertexBasePtr vertex) {
//Have 2 incoming particle and a vertex, find the possible offshell
//particles
pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
tPDSet offshells = search(vertex, inpp.first->id(), incoming,
inpp.second->id(), incoming, outgoing);
tPDSet::const_iterator it;
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isOutgoing(getParticleData(fs)) &&
vertexB->isIncoming(*it) )
final = search(vertexB, (*it)->id(), incoming, fs,
outgoing, outgoing);
//Now make diagrams
if(!final.empty())
makeDiagrams(inc, fs, final, *it, HPDiagram::sChannel,
make_pair(vertex, vertexB), make_pair(true,true));
}
}
}
void TwoToTwoProcessConstructor::
createTChannels(tPDPair inpp, long fs, tVertexBasePtr vertex) {
pair<long,long> inc = make_pair(inpp.first->id(), inpp.second->id());
//first try a with c
tPDSet offshells = search(vertex, inpp.first->id(), incoming, fs,
outgoing, outgoing);
tPDSet::const_iterator it;
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3 ) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isIncoming(inpp.second) )
final = search(vertexB, inc.second, incoming, (*it)->id(),
incoming, outgoing);
if( !final.empty() )
makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
make_pair(vertex, vertexB), make_pair(true,true));
}
}
//now try b with c
offshells = search(vertex, inpp.second->id(), incoming, fs,
outgoing, incoming);
for(it = offshells.begin(); it != offshells.end(); ++it) {
if(find(excluded_.begin(),excluded_.end(),*it)!=excluded_.end()) continue;
for(size_t iv = 0; iv < nv_; ++iv) {
tVertexBasePtr vertexB = vertices_[iv];
if( vertexB->getNpoint() != 3 ) continue;
if( !allDiagrams_ && vertexB->orderInGs() == 0 ) continue;
tPDSet final;
if( vertexB->isIncoming(inpp.first) )
final = search(vertexB, inc.first, incoming, (*it)->id(),
outgoing, outgoing);
if( !final.empty() )
makeDiagrams(inc, fs, final, *it, HPDiagram::tChannel,
make_pair(vertexB, vertex), make_pair(true, false));
}
}
}
void TwoToTwoProcessConstructor::makeFourPointDiagrams(long parta, long partb,
long partc, VBPtr vert) {
- cerr << "in make 4 " << __LINE__ << "\n";
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
getParticleData(partc));
if(loc==outgoing_.end()) return;
}
- cerr << "in make 4 " << __LINE__ << "\n";
tPDSet ext = search(vert, parta, incoming, partb,incoming, partc, outgoing);
if( ext.empty() ) return;
- cerr << "in make 4 " << __LINE__ << "\n";
IDPair in(parta, partb);
- cerr << "in make 4 " << __LINE__ << "\n";
for(tPDSet::const_iterator iter=ext.begin(); iter!=ext.end();
++iter) {
- cerr << "in make 4 " << __LINE__ << "\n";
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
*iter);
if(loc==outgoing_.end()) continue;
}
- cerr << "in make 4 " << __LINE__ << "\n";
HPDiagram nhp(in,make_pair(partc, (*iter)->id()));
nhp.vertices = make_pair(vert, vert);
nhp.channelType = HPDiagram::fourPoint;
- cerr << "in make 4 " << __LINE__ << "\n";
fixFSOrder(nhp);
- cerr << "in make 4 " << __LINE__ << "\n";
if(!checkOrder(nhp)) continue;
- cerr << "in make 4 " << __LINE__ << "\n";
if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
- cerr << "in make 4 " << __LINE__ << "\n";
}
- cerr << "in make 4 " << __LINE__ << "\n";
}
void
TwoToTwoProcessConstructor::makeDiagrams(IDPair in, long out1, const tPDSet & out2,
PDPtr inter, HPDiagram::Channel chan,
VBPair vertexpair, BPair cross) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
getParticleData(out1));
if(loc==outgoing_.end()) return;
}
for(tPDSet::const_iterator it = out2.begin(); it != out2.end(); ++it) {
if(processOption_>=1) {
PDVector::const_iterator loc = find(outgoing_.begin(),outgoing_.end(),
*it);
if(loc==outgoing_.end()) continue;
}
HPDiagram nhp( in,make_pair(out1, (*it)->id()) );
nhp.intermediate = inter;
nhp.vertices = vertexpair;
nhp.channelType = chan;
nhp.ordered = cross;
fixFSOrder(nhp);
if(!checkOrder(nhp)) continue;
if( !duplicate(nhp, processes_) ) processes_.push_back(nhp);
}
}
set<tPDPtr>
TwoToTwoProcessConstructor::search(VBPtr vertex, long part1, direction d1,
long part2, direction d2, direction d3) {
if(vertex->getNpoint() != 3) return tPDSet();
if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
vector<long> ext;
tPDSet third;
for(unsigned int ix = 0;ix < 3; ++ix) {
vector<long> pdlist = vertex->search(ix, part1);
ext.insert(ext.end(), pdlist.begin(), pdlist.end());
}
for(unsigned int ix = 0; ix < ext.size(); ix += 3) {
long id0 = ext.at(ix);
long id1 = ext.at(ix+1);
long id2 = ext.at(ix+2);
int pos;
if((id0 == part1 && id1 == part2) ||
(id0 == part2 && id1 == part1))
pos = ix + 2;
else if((id0 == part1 && id2 == part2) ||
(id0 == part2 && id2 == part1))
pos = ix + 1;
else if((id1 == part1 && id2 == part2) ||
(id1 == part2 && id2 == part1))
pos = ix;
else
pos = -1;
if(pos >= 0) {
tPDPtr p = getParticleData(ext[pos]);
if(d3 == incoming && p->CC()) p = p->CC();
third.insert(p);
}
}
return third;
}
set<tPDPtr>
TwoToTwoProcessConstructor::search(VBPtr vertex,
long part1, direction d1,
long part2, direction d2,
long part3, direction d3,
direction d4) {
if(vertex->getNpoint() != 4) return tPDSet();
if(d1 == incoming && getParticleData(part1)->CC()) part1 = -part1;
if(d2 == incoming && getParticleData(part2)->CC()) part2 = -part2;
if(d3 == incoming && getParticleData(part3)->CC()) part3 = -part3;
vector<long> ext;
tPDSet fourth;
for(unsigned int ix = 0;ix < 4; ++ix) {
vector<long> pdlist = vertex->search(ix, part1);
ext.insert(ext.end(), pdlist.begin(), pdlist.end());
}
for(unsigned int ix = 0;ix < ext.size(); ix += 4) {
long id0 = ext.at(ix); long id1 = ext.at(ix + 1);
long id2 = ext.at(ix + 2); long id3 = ext.at(ix + 3);
int pos;
if((id0 == part1 && id1 == part2 && id2 == part3) ||
(id0 == part1 && id1 == part3 && id2 == part2) ||
(id0 == part2 && id1 == part1 && id2 == part3) ||
(id0 == part2 && id1 == part3 && id2 == part1) ||
(id0 == part3 && id1 == part1 && id2 == part2) ||
(id0 == part3 && id1 == part2 && id2 == part1))
pos = ix + 3;
else if((id0 == part1 && id1 == part2 && id3 == part3) ||
(id0 == part1 && id1 == part3 && id3 == part2) ||
(id0 == part2 && id1 == part1 && id3 == part3) ||
(id0 == part2 && id1 == part3 && id3 == part1) ||
(id0 == part3 && id1 == part1 && id3 == part2) ||
(id0 == part3 && id1 == part2 && id3 == part1))
pos = ix + 2;
else if((id0 == part1 && id2 == part2 && id3 == part3) ||
(id0 == part1 && id2 == part3 && id3 == part2) ||
(id0 == part2 && id2 == part1 && id3 == part3) ||
(id0 == part2 && id2 == part3 && id3 == part1) ||
(id0 == part3 && id2 == part1 && id3 == part2) ||
(id0 == part3 && id2 == part2 && id3 == part1))
pos = ix + 1;
else if((id1 == part1 && id2 == part2 && id3 == part3) ||
(id1 == part1 && id2 == part3 && id3 == part2) ||
(id1 == part2 && id2 == part1 && id3 == part3) ||
(id1 == part2 && id2 == part3 && id3 == part1) ||
(id1 == part3 && id2 == part1 && id3 == part2) ||
(id1 == part3 && id2 == part2 && id3 == part1))
pos = ix;
else
pos = -1;
if(pos >= 0) {
tPDPtr p = getParticleData(ext[pos]);
if(d4 == incoming && p->CC())
p = p->CC();
fourth.insert(p);
}
}
return fourth;
}
void
TwoToTwoProcessConstructor::createMatrixElement(const HPDVector & process) const {
if ( process.empty() ) return;
// external particles
tcPDVector extpart(4);
extpart[0] = getParticleData(process[0].incoming.first);
extpart[1] = getParticleData(process[0].incoming.second);
extpart[2] = getParticleData(process[0].outgoing.first);
extpart[3] = getParticleData(process[0].outgoing.second);
// create the object
string objectname ("/Herwig/MatrixElements/");
string classname = MEClassname(extpart, objectname);
GeneralHardMEPtr matrixElement = dynamic_ptr_cast<GeneralHardMEPtr>
(generator()->preinitCreate(classname, objectname));
if( !matrixElement ) {
std::stringstream message;
message << "TwoToTwoProcessConstructor::createMatrixElement "
<< "- No matrix element object could be created for "
<< "the process "
<< extpart[0]->PDGName() << " " << extpart[0]->iSpin() << ","
<< extpart[1]->PDGName() << " " << extpart[1]->iSpin() << "->"
<< extpart[2]->PDGName() << " " << extpart[2]->iSpin() << ","
<< extpart[3]->PDGName() << " " << extpart[3]->iSpin()
<< ". Constructed class name: \"" << classname << "\"";
generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
return;
}
// choice for the scale
unsigned int scale;
if(scaleChoice_==0) {
// check coloured initial and final state
bool inColour = ( extpart[0]->coloured() ||
extpart[1]->coloured());
bool outColour = ( extpart[2]->coloured() ||
extpart[3]->coloured());
if(inColour&&outColour) {
bool coloured = false;
for(unsigned int ix=0;ix<process.size();++ix) {
if(process[ix].intermediate&&
process[ix].intermediate->coloured()) {
coloured = true;
break;
}
}
scale = coloured ? 1 : 0;
}
else {
scale = 0;
}
}
else {
scale = scaleChoice_-1;
}
// set the information
matrixElement->setProcessInfo(process, colourFlow(extpart),
debug(), scale, scaleFactor_);
// insert it
generator()->preinitInterface(subProcess(), "MatrixElements",
subProcess()->MEs().size(),
"insert", matrixElement->fullName());
}
string TwoToTwoProcessConstructor::MEClassname(const vector<tcPDPtr> & extpart,
string & objname) const {
string classname("Herwig::ME");
for(vector<tcPDPtr>::size_type ix = 0; ix < extpart.size(); ++ix) {
if(ix == 2) classname += "2";
if(extpart[ix]->iSpin() == PDT::Spin0) classname += "s";
else if(extpart[ix]->iSpin() == PDT::Spin1) classname += "v";
else if(extpart[ix]->iSpin() == PDT::Spin1Half) classname += "f";
else if(extpart[ix]->iSpin() == PDT::Spin2) classname += "t";
else {
std::stringstream message;
message << "MEClassname() : Encountered an unknown spin for "
<< extpart[ix]->PDGName() << " while constructing MatrixElement "
<< "classname " << extpart[ix]->iSpin();
generator()->logWarning(TwoToTwoProcessConstructorError(message.str(),Exception::warning));
}
}
objname += "ME" + extpart[0]->PDGName() + extpart[1]->PDGName() + "2"
+ extpart[2]->PDGName() + extpart[3]->PDGName();
return classname;
}