No OneTemporary
Actions

Size

35 KB

Referenced Files

None

Subscribers

None

View Options

	diff --git a/Decay/General/GeneralTwoBodyDecayer.cc b/Decay/General/GeneralTwoBodyDecayer.cc
	--- a/Decay/General/GeneralTwoBodyDecayer.cc
	+++ b/Decay/General/GeneralTwoBodyDecayer.cc
	@@ -1,367 +1,321 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the GeneralTwoBodyDecayer class.
	//

	#include "GeneralTwoBodyDecayer.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Reference.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Utilities/Exception.h"

	using namespace Herwig;

	void GeneralTwoBodyDecayer::doinit() throw(InitException) {
	DecayIntegrator::doinit();
	if(_theVertex) {
	_theVertex->init();
	}
	else
	throw InitException() << "GeneralTwoBodyDecayer::doinit() - "
	<< "Null vertex pointer!\n";

	vector<double> wgt(0);
	- PDVector inc(_theVertex->getIncoming());
	- for(unsigned int i = 0; i < inc.size(); ++i) {
	- int id = inc[i]->id();
	- if(id < 0) continue;
	- Energy m1 = inc[i]->mass();
	- PDVector decaylist(0);
	- for(unsigned int il=0;il<_thelist.size();++il) {
	- PDVector dtemp = _theVertex->search(_thelist[il],id);
	- decaylist.insert(decaylist.end(),dtemp.begin(),dtemp.end());
	+ PDVector parents = _theVertex->getIncoming();
	+ PDVector::size_type np = parents.size();
	+ PDVector extpart(3);
	+ for( tPDVector::size_type i = 0; i < np; ++i ) {
	+ tPDPtr inpart = parents[i];
	+ long pid = inpart->id();
	+ if( pid < 0 ) continue;
	+ Energy m1 = inpart->mass();
	+ PDVector decaylist;
	+ for(unsigned int il = 0; il< _thelist.size(); ++il) {
	+ PDVector temp = _theVertex->search(_thelist[il], pid);
	+ decaylist.insert(decaylist.end(),temp.begin(),temp.end());
	}
	- PDVector extpart(3);
	- for(PDVector::iterator iter=decaylist.begin();iter!=decaylist.end();) {
	- Energy m2(0.MeV),m3(0.MeV);
	- bool cc1(false),cc2(false),cc3(false);
	- if((*iter)->CC()) {cc1=true;}
	- if((*(iter+1))->CC()) {cc2=true;}
	- if((*(iter+2))->CC()) {cc3=true;}
	-
	- if((*iter)->id()==id) {
	- m2 = (*(iter+1))->mass();
	- m3 = (*(iter+2))->mass();
	- extpart[0] = *iter;
	- if(cc1) {
	- if(cc2) {extpart[1] = (*(iter+1))->CC();}
	- else {extpart[1] = *(iter+1);}
	-
	- if(cc3) {extpart[2] = (*(iter+2))->CC();}
	- else {extpart[2] = *(iter+2);}
	- }
	- else {
	- extpart[1] = *(iter+1);
	- extpart[2] = *(iter+2);
	- }
	- }
	- else if((*(iter+1))->id()==id) {
	- m2 = (*iter)->mass();
	- m3 = (*(iter+2))->mass();
	- extpart[0] = *(iter+1);
	- if(cc2) {
	- if(cc1) {extpart[1] = (*iter)->CC();}
	- else {extpart[1] = *iter;}
	-
	- if(cc3) {extpart[2] = (*(iter+2))->CC();}
	- else {extpart[2] = *(iter+2);}
	- }
	- else {
	- extpart[1] = *iter;
	- extpart[2] = *(iter+2);
	- }
	- }
	- else {
	- m2 = (*iter)->mass();
	- m3 = (*(iter+1))->mass();
	- extpart[0]=*(iter+2);
	-
	- if(cc3) {
	- if(cc1) {extpart[1] = (*iter)->CC();}
	- else {extpart[1] = *iter;}
	-
	- if(cc2) {extpart[2] = (*(iter+1))->CC();}
	- else {extpart[2] = *(iter+1);}
	- }
	- else {
	- extpart[1] = *iter;
	- extpart[2] = *(iter+1);
	- }
	-
	- }
	- if(m1 <= (m2 + m3)) {
	- decaylist.erase(iter,iter+3);
	- }
	- else {
	- _inpart.push_back(extpart[0]->id());
	- _outparta.push_back(extpart[1]->id());
	- _outpartb.push_back(extpart[2]->id());
	- DecayPhaseSpaceModePtr mode;
	- mode = new_ptr(DecayPhaseSpaceMode(extpart,this));
	- addMode(mode,_maxweight[0],wgt);
	- iter+=3;
	- }
	+ PDVector::size_type ndec = decaylist.size();
	+ for( PDVector::size_type j = 0; j < ndec; j +=3 ) {
	+ tPDPtr pa(decaylist[j]), pb(decaylist[j + 1]), pc(decaylist[j + 2]);
	+ if( pb->id() == pid ) swap(pa, pb);
	+ if( pc->id() == pid ) swap(pa, pc);
	+ //allowed on-shell decay?
	+ if( m1 <= pb->mass() + pc->mass() ) continue;
	+ //vertices are defined with all particles incoming
	+ if( pb->CC() ) pb = pb->CC();
	+ if( pc->CC() ) pc = pc->CC();
	+ //store ids so that the decayer knows what it is allowed to
	+ //decay
	+ _inpart.push_back(pid); _outparta.push_back(pb->id());
	+ _outpartb.push_back(pc->id());
	+ //create phase space mode
	+ extpart[0] = pa;
	+ extpart[1] = pb;
	+ extpart[2] = pc;
	+ addMode(new_ptr(DecayPhaseSpaceMode(extpart, this)),
	+ _maxweight[0], wgt);
	+
	}
	}
	unsigned int isize(_inpart.size()), oasize(_outparta.size()),
	obsize(_outpartb.size());
	if( isize == 0 \|\| oasize == 0 \|\| obsize == 0 )
	throw InitException()
	<< "GeneralTwoBodyDecayer::doinit() - Atleast one of the particle "
	<< "vectors has zero size, cannot continue."
	<< isize << " " << oasize << " " << obsize
	<< Exception::abortnow;

	if( isize != oasize \|\| isize != obsize )
	throw InitException()
	<< "GeneralTwoBodyDecayer::doinit() - The particle vectors have "
	<< "different sizes. " << isize << " " << oasize << " " << obsize
	<< Exception::abortnow;

	}

	int GeneralTwoBodyDecayer::modeNumber(bool & cc, tcPDPtr parent,
	const PDVector & children) const {
	long parentID = parent->id();
	long id1 = children[0]->id();
	long id2 = children[1]->id();
	int imode(-1);
	unsigned ii(0), nipart(_inpart.size());
	cc = false;
	do {
	long listpid(_inpart[ii]), listid1(_outparta[ii]),
	listid2(_outpartb[ii]);
	if( parentID == listpid &&
	((id1 == listid1 && id2 == listid2) \|\|
	(id1 == listid2 && id2 == listid1)) )
	imode = ii;
	//cc-mode
	else if(parentID == -listpid) {
	cc = true;
	if((id1 == -listid1 && id2 == -listid2) \|\|
	(id1 == -listid2 && id2 == -listid1) \|\|
	(id1 == listid1 && id2 == -listid2) \|\|
	(id1 == -listid1 && id2 == listid2) \|\|
	(id1 == listid2 && id2 == -listid1) \|\|
	(id1 == -listid2 && id2 == listid1) )
	imode = ii;
	else ++ii;
	}
	else ++ii;
	}
	while( imode < 0 && ii < nipart );
	return imode;
	}

	void GeneralTwoBodyDecayer::
	colourConnections(const Particle & parent,
	const ParticleVector & out) const {
	PDT::Colour incColour(parent.data().iColour());
	PDT::Colour outaColour(out[0]->data().iColour());
	PDT::Colour outbColour(out[1]->data().iColour());

	//incoming colour singlet
	if(incColour == PDT::Colour0) {
	// colour triplet-colourantitriplet
	if((outaColour == PDT::Colour3 && outbColour == PDT::Colour3bar) \|\|
	(outaColour == PDT::Colour3bar && outbColour == PDT::Colour3)) {
	bool ac(out[0]->id() < 0);
	out[0]->colourNeighbour(out[1],!ac);
	}
	//colour octet
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour8) {
	out[0]->colourNeighbour(out[1]);
	out[0]->antiColourNeighbour(out[1]);
	}
	// colour singlets
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour0) {
	}
	// unknown
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour singlet in "
	<< "GeneralTwoBodyDecayer::colourConnections "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	//incoming colour triplet
	else if(incColour == PDT::Colour3) {
	// colour triplet + singlet
	if(outaColour == PDT::Colour3 && outbColour == PDT::Colour0) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	//opposite order
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3) {
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	// octet + triplet
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->antiColourNeighbour(out[0]);
	}
	//opposite order
	else if(outaColour == PDT::Colour3 && outbColour == PDT::Colour8) {
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[0]->antiColourNeighbour(out[1]);
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour triplet in "
	<< "GeneralTwoBodyDecayer::colourConnections() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	// incoming colour anti triplet
	else if(incColour == PDT::Colour3bar) {
	// colour antitriplet +singlet
	if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour0) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	//opposite order
	else if(outaColour == PDT::Colour0 && outbColour == PDT::Colour3bar) {
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	//octet + antitriplet
	else if(outaColour == PDT::Colour3bar && outbColour == PDT::Colour8) {
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[0]->colourNeighbour(out[1]);
	}
	//opposite order
	else if(outaColour == PDT::Colour8 && outbColour == PDT::Colour3bar) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->colourNeighbour(out[0]);
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour antitriplet "
	<< "in GeneralTwoBodyDecayer::decay() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	//incoming colour octet
	else if(incColour == PDT::Colour8) {
	// triplet-antitriplet
	if(outaColour == PDT::Colour3&&outbColour == PDT::Colour3bar) {
	out[0]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	}
	// opposite order
	else if(outbColour == PDT::Colour3&&outaColour == PDT::Colour3bar) {
	out[0]->incomingAntiColour(const_ptr_cast<tPPtr>(&parent));
	out[1]->incomingColour(const_ptr_cast<tPPtr>(&parent));
	}
	else
	throw Exception() << "Unknown outgoing colours for decaying "
	<< "colour octet "
	<< "in GeneralTwoBodyDecayer::decay() "
	<< outaColour << " " << outbColour
	<< Exception::runerror;
	}
	else
	throw Exception() << "Unknown incoming colour in "
	<< "GeneralTwoBodyDecayer::decay() "
	<< incColour
	<< Exception::runerror;
	}

	bool GeneralTwoBodyDecayer::twoBodyMEcode(const DecayMode & dm, int & mecode,
	double & coupling) const {
	long parent = dm.parent()->id();
	ParticleMSet::const_iterator pit = dm.products().begin();
	long id1 = (*pit)->id();
	++pit;
	long id2 = (*pit)->id();
	bool order(false);
	vector<int>::size_type ix(0);
	do {
	if( parent == _inpart[ix] ) {
	long id1t(_outparta[ix]), id2t(_outpartb[ix]);
	if( id1 == id1t && id2 == id2t ) {
	order = true;
	break;
	}
	if( id1 == id2t && id2 == id1t ) {
	order = false;
	break;
	}
	}
	++ix;
	}
	while( ix < _inpart.size() );
	mecode = -1;
	coupling = 1.;
	return order;
	}


	void GeneralTwoBodyDecayer::persistentOutput(PersistentOStream & os) const {
	os << _thelist << _theVertex << _inpart << _outparta << _outpartb
	<< _maxweight;
	}

	void GeneralTwoBodyDecayer::persistentInput(PersistentIStream & is, int) {
	is >> _thelist >> _theVertex >>_inpart >>_outparta >>_outpartb
	>> _maxweight;
	}

	AbstractClassDescription<GeneralTwoBodyDecayer>
	GeneralTwoBodyDecayer::initGeneralTwoBodyDecayer;
	// Definition of the static class description member.

	void GeneralTwoBodyDecayer::Init() {

	static ClassDocumentation<GeneralTwoBodyDecayer> documentation
	("This class is designed to be a base class for all 2 body decays"
	"in a general model");

	static Reference<GeneralTwoBodyDecayer,Helicity::VertexBase> interfaceDecayVertex
	("DecayVertex",
	"Pointer to decayer vertex",
	&GeneralTwoBodyDecayer::_theVertex, false, false, true, false);

	static ParVector<GeneralTwoBodyDecayer,double> interfaceMaxWeight
	("MaxWeight",
	"Maximum weight for integration",
	&GeneralTwoBodyDecayer::_maxweight, 1.0, -1, 0, 0,
	false, false, false,&GeneralTwoBodyDecayer::setWeight, 0 ,0, 0, 0);

	static ParVector<GeneralTwoBodyDecayer,int> interfaceIncomingPart
	("IncomingPart",
	"PDG Codes for incoming particles",
	&GeneralTwoBodyDecayer::_inpart, 0, -1, 0, 0,
	false, false, false);

	static ParVector<GeneralTwoBodyDecayer,int> interfaceOutgoingPartA
	("OutgoingPartA",
	"PDG Codes for first set of outgoing particles",
	&GeneralTwoBodyDecayer::_outparta, 0, -1, 0, 0,
	false, false, false);

	static ParVector<GeneralTwoBodyDecayer,int> interfaceOutgoingPartB
	("OutgoingPartB",
	"PDG Codes for second set of outgoing particles",
	&GeneralTwoBodyDecayer::_outpartb, 0, -1, 0, 0,
	false, false, false);

	}

	double GeneralTwoBodyDecayer::brat(const DecayMode &, const Particle & p,
	double oldbrat) const {
	ParticleVector children = p.children();
	if( children.size() != 2 \|\| !p.data().widthGenerator() )
	return oldbrat;

	// partial width for this mode
	Energy scale = p.mass();
	Energy pwidth =
	partialWidth( make_pair(p.dataPtr(), scale),
	make_pair(children[0]->dataPtr(), children[0]->mass()),
	make_pair(children[1]->dataPtr(), children[1]->mass()) );
	Energy width = p.data().widthGenerator()->width(p.data(), scale);
	return pwidth/width;
	}
	diff --git a/Models/General/TwoBodyDecayConstructor.cc b/Models/General/TwoBodyDecayConstructor.cc
	--- a/Models/General/TwoBodyDecayConstructor.cc
	+++ b/Models/General/TwoBodyDecayConstructor.cc
	@@ -1,369 +1,337 @@
	// -- C++ --
	//
	// 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/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "Herwig++/Decay/General/GeneralTwoBodyDecayer.h"
	+#include "Herwig++/Models/StandardModel/StandardModel.h"

	using namespace Herwig;
	using ThePEG::Helicity::VertexBasePtr;

	+TwoBodyDecayConstructor::TwoBodyDecayConstructor():
	+ _theExistingDecayers(0),_init(true),_iteration(1),_points(1000),
	+ _info(false), _createmodes(true) {}
	+
	void TwoBodyDecayConstructor::persistentOutput(PersistentOStream & os) const {
	os << _theExistingDecayers << _init << _iteration << _points << _info;
	}

	void TwoBodyDecayConstructor::persistentInput(PersistentIStream & is, int) {
	is >>_theExistingDecayers >> _init >> _iteration >> _points >> _info;
	}

	ClassDescription<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.");

	static Switch<TwoBodyDecayConstructor,bool> interfaceInitializeDecayers
	("InitializeDecayers",
	"Initialize new decayers",
	&TwoBodyDecayConstructor::_init, true, false, false);
	static SwitchOption interfaceInitializeDecayersInitializeDecayersOn
	(interfaceInitializeDecayers,
	"On",
	"Initialize new decayers to find max weights",
	1);
	static SwitchOption interfaceInitializeDecayersoff
	(interfaceInitializeDecayers,
	"Off",
	"Use supplied weights for integration",
	0);

	static Parameter<TwoBodyDecayConstructor,int> interfaceInitIteration
	("InitIteration",
	"Number of iterations to optimise integration weights",
	&TwoBodyDecayConstructor::_iteration, 1, 0, 10,
	false, false, true);

	static Parameter<TwoBodyDecayConstructor,int> interfaceInitPoints
	("InitPoints",
	"Number of points to generate when optimising integration",
	&TwoBodyDecayConstructor::_points, 1000, 100, 100000000,
	false, false, true);

	static Switch<TwoBodyDecayConstructor,bool> interfaceOutputInfo
	("OutputInfo",
	"Whether to output information about the decayers",
	&TwoBodyDecayConstructor::_info, false, false, false);
	static SwitchOption interfaceOutputInfoOff
	(interfaceOutputInfo,
	"Off",
	"Do not output information regarding the created decayers",
	false);
	static SwitchOption interfaceOutputInfoOn
	(interfaceOutputInfo,
	"On",
	"Output information regarding the decayers",
	true);
	+
	+ static Switch<TwoBodyDecayConstructor,bool> interfaceCreateDecayModes
	+ ("CreateDecayModes",
	+ "Whether to create the ThePEG::DecayMode objects as well as the decayers",
	+ &TwoBodyDecayConstructor::_createmodes, true, false, false);
	+ static SwitchOption interfaceCreateDecayModesOn
	+ (interfaceCreateDecayModes,
	+ "On",
	+ "Create the ThePEG::DecayMode objects",
	+ true);
	+ static SwitchOption interfaceCreateDecayModesOff
	+ (interfaceCreateDecayModes,
	+ "Off",
	+ "Only create the Decayer objects",
	+ false);
	}

	-void TwoBodyDecayConstructor::DecayList(const PDVector & part) {
	- unsigned int np = part.size();
	+void TwoBodyDecayConstructor::DecayList(const PDVector & particles) {
	+ unsigned int np = particles.size();
	if( np == 0 ) return;
	- _theModel->init();
	- unsigned int nv(_theModel->numberOfVertices());
	+ tHwSMPtr model = dynamic_ptr_cast<tHwSMPtr>(generator()->standardModel());
	+ model->init();
	+ unsigned int nv(model->numberOfVertices());
	// make sure vertices are initialized
	for(unsigned int i = 0; i < nv; ++i)
	- _theModel->vertex(i)->init();
	+ model->vertex(i)->init();

	_theExistingDecayers.resize(nv,
	vector<GeneralTwoBodyDecayerPtr>(3,GeneralTwoBodyDecayerPtr()));
	- PDVector decays;
	- for(unsigned int ipart = 0; ipart < np; ++ipart) {
	+
	+ for(unsigned int ip = 0; ip < np; ++ip) {
	+ tPDPtr parent = particles[ip];
	for(unsigned int iv = 0; iv < nv; ++iv) {
	- for(unsigned int ilist = 0; ilist < 3; ++ilist) {
	- decays = createModes(part[ipart], _theModel->vertex(iv),
	- ilist, iv);
	- if(decays.size() > 0){
	- PDPtr incpart;
	- if(part[ipart]->CC())
	- incpart = part[ipart]->CC();
	- else
	- incpart = part[ipart];
	-
	- createDecayMode(incpart, decays, _theExistingDecayers[iv][ilist]);
	- }
	+ for(unsigned int il = 0; il < 3; ++il) {
	+ vector<TwoBodyDecay> decays =
	+ createModes(parent, model->vertex(iv), il, iv);
	+ if( !decays.empty() )
	+ createDecayMode(decays, _theExistingDecayers[iv][il]);
	}
	}
	}
	}

	-PDVector TwoBodyDecayConstructor::
	-createModes(tPDPtr inpart, VertexBasePtr vert,
	- unsigned int ilist, unsigned int iv) {
	+vector<TwoBodyDecay> TwoBodyDecayConstructor::
	+createModes(tPDPtr inpart, VertexBasePtr vertex,
	+ unsigned int list, unsigned int iv) {
	int id = inpart->id();
	- if(id < 0)
	- return PDVector(0);
	+ if( id < 0 \|\| !vertex->incoming(id) \|\| vertex->getNpoint() != 3 )
	+ return vector<TwoBodyDecay>();
	+ Energy m1(inpart->mass());
	+ PDVector decaylist = vertex->search(list, id);
	+ vector<TwoBodyDecay> decays;
	+ PDVector::size_type nd = decaylist.size();
	+ for( PDVector::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
	+ if( pb->CC() ) pb = pb->CC();
	+ if( pc->CC() ) pc = pc->CC();
	+ decays.push_back( make_pair(inpart, make_pair(pb, pc)) );
	+ }
	+ if( !decays.empty() )
	+ createDecayer(vertex,list,iv);

	- Energy m1(inpart->mass());
	- PDVector decaylist;
	- if(vert->getNpoint()==3 && vert->incoming(id)) {
	- decaylist = vert->search(ilist,id);
	- for(PDVector::iterator iter=decaylist.begin();iter!=decaylist.end();) {
	- Energy m2,m3;
	- bool cc1(false),cc2(false),cc3(false);
	- if((*iter)->CC()) cc1=true;
	- if((*(iter+1))->CC()) cc2=true;
	- if((*(iter+2))->CC()) cc3=true;
	-
	- if((*iter)->id()==id) {
	- m2 = (*(iter+1))->mass();
	- m3 = (*(iter+2))->mass();
	- if(cc1) {
	- if(cc2)
	- (iter+1) = ((iter+1))->CC();
	- if(cc3)
	- (iter+2) = ((iter+2))->CC();
	- }
	- }
	- else if((*(iter+1))->id()==id) {
	- m2 = (*iter)->mass();
	- m3 = (*(iter+2))->mass();
	- if(cc2) {
	- if(cc1)
	- iter = (iter)->CC();
	- if(cc3)
	- (iter+2) = ((iter+2))->CC();
	- }
	- }
	- else {
	- m2 = (*iter)->mass();
	- m3 = (*(iter+1))->mass();
	- if(cc3) {
	- if(cc1)
	- iter = (iter)->CC();
	- if(cc2)
	- (iter+1) = ((iter+1))->CC();
	- }
	- }
	-
	- if(m1 <= (m2 + m3))
	- decaylist.erase(iter,iter+3);
	- else
	- iter+=3;
	- }
	-
	- if(decaylist.size() > 0)
	- createDecayer(vert,ilist,iv);
	- }
	-
	- return decaylist;
	+ return decays;
	}

	-void TwoBodyDecayConstructor::createDecayer(VertexBasePtr vert,
	+void TwoBodyDecayConstructor::createDecayer(VertexBasePtr vertex,
	unsigned int icol,
	unsigned int ivert) {
	if( _theExistingDecayers[ivert][icol] ) return;
	string name;
	- switch(vert->getName()) {
	+ switch(vertex->getName()) {
	case FFV : {
	if( icol == 0 \|\| icol == 1)
	name = "FFVDecayer";
	else
	name = "VFFDecayer";
	}
	break;
	case FFS : {
	if( icol == 0 \|\| icol == 1)
	name = "FFSDecayer";
	else
	name = "SFFDecayer";
	}
	break;
	case VVS : {
	if( icol == 0 \|\| icol == 1)
	name = "VVSDecayer";
	else
	name = "SVVDecayer";
	}
	break;
	case GeneralSVV : {
	if(icol == 0)
	name = "SVVLoopDecayer";
	}
	break;
	case VSS : {
	if(icol == 0)
	name = "VSSDecayer";
	else
	name = "SSVDecayer";
	}
	break;
	case VVT : {
	if(icol == 2)
	name = "TVVDecayer";
	}
	break;
	case FFT : {
	if(icol == 2)
	name = "TFFDecayer";
	}
	break;
	case SST : {
	if(icol == 2)
	name = "TSSDecayer";
	}
	break;
	case SSS : {
	name = "SSSDecayer";
	}
	break;
	case VVV : name = "VVVDecayer";
	break;
	default : throw NBodyDecayConstructorError()
	- << "Error: Cannot assign " << vert->fullName() << " to a decayer. "
	+ << "Error: Cannot assign " << vertex->fullName() << " to a decayer. "
	<< "Looking in column " << icol;
	}
	ostringstream fullname;
	fullname << "/Defaults/Decays/" << name << "_"
	<< ivert << "_" << icol;
	string classname = "Herwig::" + name;

	GeneralTwoBodyDecayerPtr decayer;
	decayer = dynamic_ptr_cast<GeneralTwoBodyDecayerPtr>
	(generator()->preinitCreate(classname,fullname.str()));
	string msg = generator()->preinitInterface(decayer, "DecayVertex",
	- "set", vert->fullName());
	+ "set", vertex->fullName());
	if(msg.find("Error:") != string::npos)
	throw NBodyDecayConstructorError()
	<< "TwoBodyDecayConstructor::createDecayer - An error occurred while "
	<< "setting the vertex for " << decayer->fullName()
	<< " - " << msg
	<< Exception::abortnow;
	decayer->init();
	setDecayerInterfaces(fullname.str());
	_theExistingDecayers[ivert][icol] = decayer;
	}

	void TwoBodyDecayConstructor::
	-createDecayMode(tPDPtr inpart, const PDVector & decays,
	+createDecayMode(const vector<TwoBodyDecay> & decays,
	GeneralTwoBodyDecayerPtr decayer) {
	if(!decayer)
	throw NBodyDecayConstructorError()
	<< "TwoBodyDecayConstructor::createDecayMode - The decayer "
	<< "pointer is null!\n"
	<< Exception::runerror;
	- PDVector children(2);
	- if(inpart->CC())
	- inpart = (inpart->CC());
	+ tPDPtr inpart = decays[0].first;
	inpart->stable(false);
	tEGPtr eg = generator();
	- for(unsigned int ix = 0; ix < decays.size(); ix += 3) {
	- if(decays[ix]->id() == inpart->id()) {
	- children[0] = decays[ix+1];
	- children[1] = decays[ix+2];
	- }
	- else if(decays[ix+1]->id() == inpart->id()) {
	- children[0] = decays[ix];
	- children[1] = decays[ix+2];
	- }
	- else {
	- children[0] = decays[ix];
	- children[1] = decays[ix+1];
	- }
	- string tag = inpart->PDGName() + "->" + children[0]->PDGName() +
	- "," + children[1]->PDGName() + ";";
	+ for(unsigned int ix = 0; ix < decays.size(); ++ix ) {
	+ tPDPtr pb(decays[ix].second.first), pc(decays[ix].second.second);
	+ string tag = inpart->PDGName() + "->" + pb->PDGName() +
	+ "," + pc->PDGName() + ";";
	//now create DecayMode objects that do not already exist
	tDMPtr dm = eg->findDecayMode(tag);
	if ( !dm ) {
	- tag = inpart->PDGName() + "->" + children[1]->PDGName() +
	- "," + children[0]->PDGName() + ";";
	+ tag = inpart->PDGName() + "->" + pc->PDGName() +
	+ "," + pb->PDGName() + ";";
	dm = eg->findDecayMode(tag);
	}
	- if( !dm ) {
	+ if( !dm && _createmodes ) {
	tDMPtr ndm = eg->preinitCreateDecayMode(tag);
	if(ndm) {
	eg->preinitInterface(ndm, "Decayer", "set",
	decayer->fullName());
	eg->preinitInterface(ndm, "OnOff", "set", "1");
	Energy width =
	decayer->partialWidth(make_pair(inpart,inpart->mass()),
	- make_pair(children[0],children[0]->mass()) ,
	- make_pair(children[1], children[1]->mass()));
	+ make_pair(pb,pb->mass()) ,
	+ make_pair(pc,pc->mass()));
	setBranchingRatio(ndm, width);
	}
	else
	throw NBodyDecayConstructorError()
	<< "TwoBodyDecayConstructor::createDecayMode - Needed to create "
	<< "new decaymode but one could not be created for the tag "
	<< tag << Exception::warning;
	}
	- else {
	+ else if( dm ) {
	if((dm->decayer()->fullName()).find("Mambo") != string::npos)
	eg->preinitInterface(dm, "Decayer", "set",
	decayer->fullName());
	}
	+ else {}
	}
	//update CC mode if it exists
	if( inpart->CC() )
	inpart->CC()->synchronize();
	}

	void TwoBodyDecayConstructor::setBranchingRatio(tDMPtr dm, Energy pwidth) {
	//Need width and branching ratios for all currently created decay modes
	PDPtr parent = const_ptr_cast<PDPtr>(dm->parent());
	Selector<tDMPtr> modes = parent->decaySelector();
	Energy currentwidth(0.*MeV);
	if( !modes.empty() ) currentwidth = parent->width();
	Energy newWidth = currentwidth + pwidth;
	parent->width(newWidth);
	//need to reweight current branching fractions if there are any
	for(Selector<tDMPtr>::const_iterator dit = modes.begin();
	dit != modes.end(); ++dit) {
	double newbrat = ((dit).second->brat())currentwidth/newWidth;
	ostringstream brf;
	brf << newbrat;
	generator()->preinitInterface((*dit).second, "BranchingRatio",
	"set", brf.str());
	}
	//set brat for current mode
	double brat = pwidth/newWidth;
	ostringstream br;
	br << brat;
	generator()->preinitInterface(dm, "BranchingRatio",
	"set", br.str());
	parent->touch();
	parent->update();
	parent->reset();
	}

	void TwoBodyDecayConstructor::setDecayerInterfaces(string fullname) const {
	if( _init ) {
	ostringstream value;
	value << _init;
	generator()->preinitInterface(fullname, "Initialize", "set",
	value.str());
	value.str("");
	value << _iteration;
	generator()->preinitInterface(fullname, "Iteration", "set",
	value.str());
	value.str("");
	value << _points;
	generator()->preinitInterface(fullname, "Points", "set",
	value.str());
	}
	string outputmodes;
	if( _info ) outputmodes = string("Output");
	else outputmodes = string("NoOutput");
	generator()->preinitInterface(fullname, "OutputModes", "set",
	outputmodes);
	}
	diff --git a/Models/General/TwoBodyDecayConstructor.h b/Models/General/TwoBodyDecayConstructor.h
	--- a/Models/General/TwoBodyDecayConstructor.h
	+++ b/Models/General/TwoBodyDecayConstructor.h
	@@ -1,236 +1,225 @@
	// -- C++ --
	#ifndef HERWIG_TwoBodyDecayConstructor_H
	#define HERWIG_TwoBodyDecayConstructor_H
	//
	// This is the declaration of the TwoBodyDecayConstructor class.
	//

	#include "NBodyDecayConstructorBase.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	-#include "Herwig++/Decay/DecayIntegrator.h"
	-#include "Herwig++/Decay/DecayPhaseSpaceMode.h"
	#include "Herwig++/Decay/General/GeneralTwoBodyDecayer.fh"
	-#include "Herwig++/Models/StandardModel/StandardModel.h"
	#include "TwoBodyDecayConstructor.fh"

	namespace Herwig {
	using namespace ThePEG;
	using Helicity::VertexType;
	using Helicity::VertexBasePtr;

	/**
	+ * A two body decay mode
	+ */
	+typedef pair<tPDPtr, tPDPair> TwoBodyDecay;
	+
	+/**
	* The TwoBodyDecayConstructor class inherits from the dummy base class
	* NBodyDecayConstructorBase and implements the necessary functions in
	* order to create the 2 body decaymodes for a given set of vertices
	* stored in a Model class.
	*
	* @see \ref TwoBodyDecayConstructorInterfaces "The interfaces"
	* defined for TwoBodyDecayConstructor.
	* @see NBodyDecayConstructor
	**/
	class TwoBodyDecayConstructor: public NBodyDecayConstructorBase {
	-
	+
	public:

	/**
	* The default constructor.
	*/
	- inline TwoBodyDecayConstructor();
	+ TwoBodyDecayConstructor();

	/**
	* Function used to determine allowed decaymodes
	*@param part vector of ParticleData pointers containing particles in model
	*/
	virtual void DecayList(const PDVector & part);

	public:

	/** @name Functions used by the persistent I/O system. */
	//@{
	/**
	* Function used to write out object persistently.
	* @param os the persistent output stream written to.
	*/
	void persistentOutput(PersistentOStream & os) const;

	/**
	* Function used to read in object persistently.
	* @param is the persistent input stream read from.
	* @param version the version number of the object when written.
	*/
	void persistentInput(PersistentIStream & is, int version);
	//@}

	/**
	* The standard Init function used to initialize the interfaces.
	* Called exactly once for each class by the class description system
	* before the main function starts or
	* when this class is dynamically loaded.
	*/
	static void Init();

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	inline virtual IBPtr clone() const;

	/** Make a clone of this object, possibly modifying the cloned object
	* to make it sane.
	* @return a pointer to the new object.
	*/
	inline virtual IBPtr fullclone() const;
	//@}

	-protected:
	-
	- /** @name Standard Interfaced functions. */
	- //@{
	- /**
	- * Initialize this object after the setup phase before saving an
	- * EventGenerator to disk.
	- * @throws InitException if object could not be initialized properly.
	- */
	- inline virtual void doinit() throw(InitException);
	- //@}
	-
	private:

	/**
	* The static object used to initialize the description of this class.
	* Indicates that this is a concrete class with persistent data.
	*/
	static ClassDescription<TwoBodyDecayConstructor> initTwoBodyDecayConstructor;

	/**
	* The assignment operator is private and must never be called.
	* In fact, it should not even be implemented.
	*/
	TwoBodyDecayConstructor & operator=(const TwoBodyDecayConstructor &);

	private:

	/** @name Functions to create decayers and decaymodes. */
	//@{
	/**
	* Function to create decays
	* @param inpart Incoming particle
	* @param vert The vertex to create decays for
	* @param ilist Which list to search
	* @param iv Row number in _theExistingDecayers member
	- * @return vector of ParticleData ptrs
	+ * @return A vector a decay modes
	*/
	- PDVector createModes(tPDPtr inpart, VertexBasePtr vert,
	- unsigned int ilist,
	- unsigned int iv);
	+ vector<TwoBodyDecay> createModes(tPDPtr inpart, VertexBasePtr vert,
	+ unsigned int ilist,
	+ unsigned int iv);

	/**
	* Function to create decayer for specific vertex
	* @param vert Pointer to vertex
	* @param icol Integer referring to the colmun in _theExistingDecayers
	* @param ivert Integer referring to the row in _theExistingDecayers
	* member variable
	*/
	void createDecayer(VertexBasePtr vert, unsigned int icol,
	unsigned int ivert);

	/**
	* Create decay mode(s) from given part and decay modes
	- * @param inpart pointer to incoming particle
	- * @param decays list of allowed interactions
	+ * @param decays The vector of decay modes
	* @param decayer The decayer responsible for this decay
	*/
	- void createDecayMode(tPDPtr inpart,
	- const PDVector & decays,
	+ void createDecayMode(const vector<TwoBodyDecay> & decays,
	GeneralTwoBodyDecayerPtr decayer);

	/**
	* Set the branching ratio of this mode. This requires
	* calculating a new width for the decaying particle and reweighting
	* the current branching fractions.
	* @param dm The decaymode for which to set the branching ratio
	* @param pwidth The calculated width of the mode
	*/
	void setBranchingRatio(tDMPtr dm, Energy pwidth);

	/**
	* Set the interfaces of the decayers depending on the flags stored.
	* @param name Fullname of the decayer in the EventGenerator
	* including the path
	*/
	void setDecayerInterfaces(string name) const;
	//@}

	private:

	/**
	* Existing decayers
	*/
	vector<vector<GeneralTwoBodyDecayerPtr> > _theExistingDecayers;

	/**
	- * Model Pointer
	- */
	- Ptr<Herwig::StandardModel>::pointer _theModel;
	-
	- /**
	* Whether to initialize decayers or not
	*/
	bool _init;

	/**
	* Number of iterations if initializing (default 1)
	*/
	int _iteration;

	/**
	* Number of points to do in initialization
	*/
	int _points;

	/**
	* Whether to output information on the decayers
	*/
	bool _info;
	+
	+
	+ /**
	+ * Whether to create the DecayModes as well as the Decayer objects
	+ */
	+ bool _createmodes;
	};

	}

	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/** This template specialization informs ThePEG about the
	* base classes of TwoBodyDecayConstructor. */
	template <>
	struct BaseClassTrait<Herwig::TwoBodyDecayConstructor,1> {
	/** Typedef of the first base class of TwoBodyDecayConstructor. */
	typedef Herwig::NBodyDecayConstructorBase NthBase;
	};

	/** This template specialization informs ThePEG about the name of
	* the TwoBodyDecayConstructor class and the shared object where it is defined. */
	template <>
	struct ClassTraits<Herwig::TwoBodyDecayConstructor>
	: public ClassTraitsBase<Herwig::TwoBodyDecayConstructor> {
	/** Return a platform-independent class name */
	static string className() { return "Herwig::TwoBodyDecayConstructor"; }
	/** Return the name of the shared library be loaded to get
	* access to the TwoBodyDecayConstructor class and every other class it uses
	* (except the base class). */
	static string library() { return "libHwModelGenerator.so"; }
	};

	/** @endcond */

	}

	#include "TwoBodyDecayConstructor.icc"

	#endif /* HERWIG_TwoBodyDecayConstructor_H */
	diff --git a/Models/General/TwoBodyDecayConstructor.icc b/Models/General/TwoBodyDecayConstructor.icc
	--- a/Models/General/TwoBodyDecayConstructor.icc
	+++ b/Models/General/TwoBodyDecayConstructor.icc
	@@ -1,27 +1,17 @@
	// -- C++ --
	//
	// This is the implementation of the inlined member functions of
	// the TwoBodyDecayConstructor class.
	//

	namespace Herwig {
	-
	-inline TwoBodyDecayConstructor::TwoBodyDecayConstructor():
	- _theExistingDecayers(0),_init(true),_iteration(1),_points(1000),
	- _info(false) {}

	inline IBPtr TwoBodyDecayConstructor::clone() const {
	return new_ptr(*this);
	}

	inline IBPtr TwoBodyDecayConstructor::fullclone() const {
	return new_ptr(*this);
	}

	-inline void TwoBodyDecayConstructor::doinit() throw(InitException) {
	- NBodyDecayConstructorBase::doinit();
	- _theModel = dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>
	- (generator()->standardModel());
	}
	-
	-}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Sep 30, 6:16 AM (1 h, 23 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 6553914
Default Alt Text: (35 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions