No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

	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();
	}

File Metadata

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

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions