No OneTemporary
Actions

Size

17 KB

Subscribers

None

View Options

	diff --git a/Models/General/WeakCurrentDecayConstructor.cc b/Models/General/WeakCurrentDecayConstructor.cc
	--- a/Models/General/WeakCurrentDecayConstructor.cc
	+++ b/Models/General/WeakCurrentDecayConstructor.cc
	@@ -1,296 +1,295 @@
	// -- C++ --
	//
	// WeakCurrentDecayConstructor.cc is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 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 WeakCurrentDecayConstructor class.
	//

	#include "WeakCurrentDecayConstructor.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Interface/RefVector.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/PDT/DecayMode.h"
	#include "ThePEG/Helicity/Vertex/AbstractFFVVertex.fh"
	#include "DecayConstructor.h"

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

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

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

	void WeakCurrentDecayConstructor::doinit() {
	NBodyDecayConstructorBase::doinit();
	- _theModel = dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>
	- (generator()->standardModel());
	+ model_ = dynamic_ptr_cast<Ptr<Herwig::StandardModel>::pointer>(generator()->standardModel());
	unsigned int isize=decayTags_.size();
	if(isize!=_norm .size()\|\|isize!=_current.size())
	throw InitException() << "Invalid sizes for the decay mode vectors in "
	<< " WeakCurrentDecayConstructor "
	<< decayTags_.size() << " " << _norm.size() << " "
	<< _current.size() << Exception::runerror;
	// get the particles from the tags
	for(unsigned int ix=0;ix<decayTags_.size();++ix) {
	_current[ix]->init();
	particles_.push_back(vector<tPDPtr>());
	string tag=decayTags_[ix];
	do {
	string::size_type next = min(tag.find(','), tag.find(';'));
	particles_.back().push_back(generator()->findParticle(tag.substr(0,next)));
	if(!particles_.back().back())
	throw Exception() << "Failed to find particle " << tag.substr(0,next)
	<< " in DecayMode " << decayTags_[ix]
	<< " in WeakCurrentDecayConstructor::doinit()"
	<< Exception::runerror;
	if(tag[next]==';') break;
	tag = tag.substr(next+1);
	}
	while(true);
	}
	}

	void WeakCurrentDecayConstructor::persistentOutput(PersistentOStream & os) const {
	- os << ounit(_masscut,GeV) << decayTags_ << particles_ << _norm << _current;
	+ os << ounit(massCut_,GeV) << decayTags_ << particles_ << _norm << _current;
	}

	void WeakCurrentDecayConstructor::persistentInput(PersistentIStream & is, int) {
	- is >> iunit(_masscut,GeV) >> decayTags_ >> particles_ >> _norm >> _current;
	+ is >> iunit(massCut_,GeV) >> decayTags_ >> particles_ >> _norm >> _current;
	}

	// The following static variable is needed for the type
	// description system in ThePEG.
	DescribeClass<WeakCurrentDecayConstructor,NBodyDecayConstructorBase>
	describeHerwigWeakCurrentDecayConstructor("Herwig::WeakCurrentDecayConstructor", "Herwig.so");

	void WeakCurrentDecayConstructor::Init() {

	static ClassDocumentation<WeakCurrentDecayConstructor> documentation
	("The WeakCurrentDecayConstructor class implemets the decay of BSM particles "
	"to low mass hadronic states using the Weak current");

	static ParVector<WeakCurrentDecayConstructor,string> interfaceDecayModes
	("DecayModes",
	"The decays of the weak current",
	&WeakCurrentDecayConstructor::decayTags_, -1, "", "", "",
	false, false, Interface::nolimits);

	static ParVector<WeakCurrentDecayConstructor,double> interfaceNormalisation
	("Normalisation",
	"The normalisation of the different modes",
	&WeakCurrentDecayConstructor::_norm, -1, 1.0, 0.0, 10.0,
	false, false, Interface::limited);

	static RefVector<WeakCurrentDecayConstructor,WeakCurrent> interfaceCurrent
	("Current",
	"The current for the decay mode",
	&WeakCurrentDecayConstructor::_current, -1, false, false, true, false, false);

	static Parameter<WeakCurrentDecayConstructor,Energy> interfaceMassCut
	("MassCut",
	"The maximum mass difference for the decay",
	- &WeakCurrentDecayConstructor::_masscut, GeV, 5.0GeV, 1.0GeV, 10.0*GeV,
	+ &WeakCurrentDecayConstructor::massCut_, GeV, 5.0GeV, 1.0GeV, 10.0*GeV,
	false, false, Interface::limited);

	}

	void WeakCurrentDecayConstructor::DecayList(const set<PDPtr> & part) {
	if( part.empty() ) return;
	- unsigned int nv(_theModel->numberOfVertices());
	+ unsigned int nv(model_->numberOfVertices());
	for(set<PDPtr>::const_iterator ip=part.begin();ip!=part.end();++ip) {
	for(unsigned int iv = 0; iv < nv; ++iv) {
	for(unsigned int ilist = 0; ilist < 3; ++ilist) {
	vector<TwoBodyDecay> decays =
	- createModes(*ip, _theModel->vertex(iv),ilist);
	+ createModes(*ip, model_->vertex(iv),ilist);
	if(!decays.empty()) createDecayMode(decays);
	}
	}
	}
	}

	vector<TwoBodyDecay> WeakCurrentDecayConstructor::createModes(const PDPtr inpart,
	const VertexBasePtr vert,
	unsigned int ilist) {
	int id = inpart->id();
	if( !vert->isIncoming(inpart) \|\| vert->getNpoint() != 3 )
	return vector<TwoBodyDecay>();
	Energy m1(inpart->mass());
	vector<tPDPtr> decaylist;
	decaylist = vert->search(ilist,inpart);
	tPDVector::size_type nd = decaylist.size();
	vector<TwoBodyDecay> decays;
	for( tPDVector::size_type i = 0; i < nd; i += 3 ) {
	tPDPtr pa(decaylist[i]), pb(decaylist.at(i + 1)),
	pc(decaylist.at(i + 2));
	if( pb->id() == id ) swap(pa, pb);
	if( pc->id() == id ) swap(pa, pc);
	//One of the products must be a W
	Energy mp(ZERO);
	if( abs(pb->id()) == ParticleID::Wplus )
	mp = pc->mass();
	else if( abs(pc->id()) == ParticleID::Wplus )
	mp = pb->mass();
	else
	continue;
	//allowed on-shell decay and passes mass cut
	if( m1 >= pb->mass() + pc->mass() ) continue;
	if( m1 < mp ) continue;
	- if( m1 - mp >= _masscut ) continue;
	+ if( m1 - mp >= massCut_ ) continue;
	//vertices are defined with all particles incoming
	if( pb->CC() ) pb = pb->CC();
	if( pc->CC() ) pc = pc->CC();
	decays.push_back( TwoBodyDecay(inpart,pb, pc, vert) );
	if(abs(decays.back().children_.second->id())!=ParticleID::Wplus)
	swap(decays.back().children_.first,decays.back().children_.second);
	assert(abs(decays.back().children_.second->id())==ParticleID::Wplus);
	}
	return decays;
	}

	GeneralCurrentDecayerPtr WeakCurrentDecayConstructor::createDecayer(PDPtr in, PDPtr out1,
	vector<tPDPtr> outCurrent,
	VertexBasePtr vertex,
	WeakCurrentPtr current) {
	string name;
	using namespace ThePEG::Helicity::VertexType;
	switch(vertex->getName()) {
	case FFV :
	name = "FFVCurrentDecayer";
	break;
	default :
	ostringstream message;
	message << "Invalid vertex for decays of " << in->PDGName() << " -> " << out1->PDGName()
	<< " via weak current " << vertex->fullName() << "\n";
	generator()->logWarning(NBodyDecayConstructorError(message.str(),
	Exception::warning));
	return GeneralCurrentDecayerPtr();
	}
	ostringstream fullname;
	fullname << "/Herwig/Decays/" << name << "_" << in->PDGName() << "_"
	<< out1->PDGName();
	for(unsigned int ix=0;ix<outCurrent.size();++ix)
	fullname << "_" << outCurrent[ix]->PDGName();
	string classname = "Herwig::" + name;
	GeneralCurrentDecayerPtr decayer = dynamic_ptr_cast<GeneralCurrentDecayerPtr>
	(generator()->preinitCreate(classname,fullname.str()));
	- decayer->setDecayInfo(in,out1,outCurrent,vertex,current,_masscut);
	+ decayer->setDecayInfo(in,out1,outCurrent,vertex,current,massCut_);
	// set decayer options from base class
	setDecayerInterfaces(fullname.str());
	// initialize the decayer
	decayer->init();
	// return the decayer
	return decayer;
	}

	void WeakCurrentDecayConstructor::
	createDecayMode(vector<TwoBodyDecay> & decays) {
	assert(!decays.empty());
	for(unsigned int ix = 0; ix < decays.size(); ++ix) {
	PDVector particles(3);
	particles[0] = decays[ix].parent_;
	particles[1] = decays[ix].children_.first ;
	bool Wplus=decays[ix].children_.second->id()==ParticleID::Wplus;
	for(unsigned int iy=0;iy<_current.size();++iy) {
	particles.resize(2);
	vector<tPDPtr> wprod=particles_[iy];
	int icharge=0;
	Energy msum = particles[0]->mass()-particles[1]->mass();
	for(unsigned int iz=0;iz<wprod.size();++iz) {
	icharge += wprod[iz]->iCharge();
	msum -=wprod[iz]->mass();
	}
	if(msum<=ZERO) continue;
	bool cc = (Wplus&&icharge==-3)\|\|(!Wplus&&icharge==3);
	OrderedParticles outgoing;
	outgoing.insert(particles[1]);
	for(unsigned int iz=0;iz<wprod.size();++iz) {
	if(cc&&wprod[iz]->CC()) wprod[iz]=wprod[iz]->CC();
	outgoing.insert(wprod[iz]);
	}
	// check outgoing particles initialised
	for(unsigned int iz=0;iz<wprod.size();++iz) wprod[iz]->init();
	// create the tag for the decay mode
	string tag = particles[0]->PDGName() + "->";
	OrderedParticles::const_iterator it = outgoing.begin();
	do {
	tag += (**it).name();
	++it;
	if(it!=outgoing.end()) tag +=",";
	else tag +=";";
	}
	while(it!=outgoing.end());
	// find the decay mode
	tDMPtr dm= generator()->findDecayMode(tag);
	if( !dm && createDecayModes() ) {
	// create the decayer
	GeneralCurrentDecayerPtr decayer = createDecayer(particles[0],particles[1],
	wprod,decays[ix].vertex_,
	_current[iy]);
	if(!decayer) continue;
	// calculate the width
	Energy pWidth = _norm[iy]*decayer->partialWidth(particles[0],particles[1],wprod);
	if(pWidth<=ZERO) {
	generator()->preinitInterface(decayer->fullName(),
	"Initialize", "set","0");
	continue;
	}
	tDMPtr ndm = generator()->preinitCreateDecayMode(tag);
	if(!ndm) throw NBodyDecayConstructorError()
	<< "WeakCurrentDecayConstructor::createDecayMode - Needed to create "
	<< "new decaymode but one could not be created for the tag "
	<< tag
	<< Exception::warning;
	generator()->preinitInterface(ndm, "Decayer", "set",
	decayer->fullName());
	generator()->preinitInterface(ndm, "Active", "set", "Yes");
	setBranchingRatio(ndm, pWidth);
	particles[0]->stable(false);
	if(ndm->brat()<decayConstructor()->minimumBR()) {
	generator()->preinitInterface(decayer->fullName(),
	"Initialize", "set","0");
	}
	}
	else if (dm) {
	// create the decayer
	GeneralCurrentDecayerPtr decayer = createDecayer(particles[0],particles[1],
	wprod,decays[ix].vertex_,
	_current[iy]);
	if(!decayer) continue;
	generator()->preinitInterface(dm, "Decayer", "set", decayer->fullName());
	particles[0]->stable(false);
	if(createDecayModes()) {
	// calculate the width
	Energy pWidth = _norm[iy]*decayer->partialWidth(particles[0],particles[1],wprod);
	if(pWidth<=ZERO) {
	generator()->preinitInterface(decayer->fullName(),
	"Initialize", "set","0");
	continue;
	}
	generator()->preinitInterface(dm, "Active", "set", "Yes");
	particles[0]->width(particles[0]->width()*(1.-dm->brat()));
	setBranchingRatio(dm, pWidth);
	}
	if(dm->brat()<decayConstructor()->minimumBR()) {
	generator()->preinitInterface(decayer->fullName(),
	"Initialize", "set","0");
	}
	}
	}
	}
	}
	diff --git a/Models/General/WeakCurrentDecayConstructor.h b/Models/General/WeakCurrentDecayConstructor.h
	--- a/Models/General/WeakCurrentDecayConstructor.h
	+++ b/Models/General/WeakCurrentDecayConstructor.h
	@@ -1,194 +1,194 @@
	// -- C++ --
	//
	// WeakCurrentDecayConstructor.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2019 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.
	//
	#ifndef HERWIG_WeakCurrentDecayConstructor_H
	#define HERWIG_WeakCurrentDecayConstructor_H
	//
	// This is the declaration of the WeakCurrentDecayConstructor class.
	//

	#include "NBodyDecayConstructorBase.h"
	#include "ThePEG/Helicity/Vertex/VertexBase.h"
	#include "Herwig/Decay/General/GeneralCurrentDecayer.fh"
	#include "Herwig/Models/StandardModel/StandardModel.h"
	#include "Herwig/Decay/WeakCurrents/WeakCurrent.h"
	#include "Herwig/Decay/General/GeneralCurrentDecayer.h"
	#include "TwoBodyDecay.h"

	namespace Herwig {

	using namespace ThePEG;

	/**
	* Here is the documentation of the WeakCurrentDecayConstructor class.
	*
	* @see \ref WeakCurrentDecayConstructorInterfaces "The interfaces"
	* defined for WeakCurrentDecayConstructor.
	*/
	class WeakCurrentDecayConstructor: public NBodyDecayConstructorBase {

	public:

	/**
	* The default constructor.
	*/
	- WeakCurrentDecayConstructor() : _masscut(5.*GeV) {}
	+ WeakCurrentDecayConstructor() : massCut_(5.*GeV) {}

	/**
	* Function used to determine allowed decaymodes, to be implemented
	* in derived class.
	*@param part vector of ParticleData pointers containing particles in model
	*/
	virtual void DecayList(const set<PDPtr> & part);

	/**
	* Number of outgoing lines. Required for correct ordering (do this one last)
	*/
	virtual unsigned int numBodies() const { return 1000; }

	/**
	* Cut off
	*/
	- Energy massCut() const { return _masscut;}
	+ Energy massCut() const { return massCut_;}

	public:

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

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

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

	protected:

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

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

	protected:

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

	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
	*/
	vector<TwoBodyDecay>
	createModes(const PDPtr inpart,const VertexBasePtr vert,
	unsigned int ilist);

	/**
	* 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
	*/
	GeneralCurrentDecayerPtr createDecayer(PDPtr in, PDPtr out1,
	vector<tPDPtr> outCurrent,
	VertexBasePtr vertex,
	WeakCurrentPtr current);

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

	private:

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

	private:

	/**
	* Model Pointer
	*/
	- Ptr<Herwig::StandardModel>::pointer _theModel;
	+ Ptr<Herwig::StandardModel>::pointer model_;

	/**
	* Cut-off on the mass difference
	*/
	- Energy _masscut;
	+ Energy massCut_;

	/**
	* Tags for the modes
	*/
	vector<string> decayTags_;

	/**
	* Particles for the mode
	*/
	vector<vector<tPDPtr> > particles_;

	/**
	* Normalisation
	*/
	vector<double> _norm;

	/**
	* The current for the mode
	*/
	vector<WeakCurrentPtr> _current;
	};

	}

	#endif /* HERWIG_WeakCurrentDecayConstructor_H */

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:19 PM (1 d, 16 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3795742
Default Alt Text: (17 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions