No OneTemporary
Actions

Size

33 KB

Subscribers

None

View Options

	diff --git a/Models/ADD/GravitonMassGenerator.cc b/Models/ADD/GravitonMassGenerator.cc
	--- a/Models/ADD/GravitonMassGenerator.cc
	+++ b/Models/ADD/GravitonMassGenerator.cc
	@@ -1,91 +1,92 @@
	// -- C++ --
	//
	// This is the implementation of the non-inlined, non-templated member
	// functions of the GravitonMassGenerator class.
	//

	#include "GravitonMassGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	+#include "Herwig/PDT/GenericWidthGenerator.h"
	#include "ADDModel.h"

	using namespace Herwig;

	GravitonMassGenerator::GravitonMassGenerator()
	: prefactor_(0.), delta_(2), md_(1000.*GeV), mMin_(MeV)
	{}

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

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

	void GravitonMassGenerator::persistentOutput(PersistentOStream & os) const {
	os << prefactor_ << delta_ << ounit(md_,GeV) << ounit(mMin_,GeV);
	}

	void GravitonMassGenerator::persistentInput(PersistentIStream & is, int) {
	is >> prefactor_ >> delta_ >> iunit(md_,GeV) >> iunit(mMin_,GeV);
	}

	ClassDescription<GravitonMassGenerator>
	GravitonMassGenerator::initGravitonMassGenerator;
	// Definition of the static class description member.

	void GravitonMassGenerator::Init() {

	static ClassDocumentation<GravitonMassGenerator> documentation
	("The GravitonMassGenerator class generates the mass for external gravitions "
	"in the ADD model.");

	static Parameter<GravitonMassGenerator,Energy> interfaceMinimumMass
	("MinimumMass",
	"Minimum gravition mass to avoid numerical problems",
	&GravitonMassGenerator::mMin_, GeV, MeV, MeV, GeV,
	false, false, Interface::limited);

	}

	void GravitonMassGenerator::doinit() {
	GenericMassGenerator::doinit();
	tcHwADDPtr hwADD = dynamic_ptr_cast<tcHwADDPtr>(generator()->standardModel());
	if(!hwADD)
	throw Exception() << "Must have ADDModel in GravitonMassGenerator::doinit()"
	<< Exception::runerror;
	delta_ = hwADD->delta();
	md_ = hwADD->MD();
	// calculate the prefactor
	prefactor_ = sqr(hwADD->MPlanckBar()/md_);
	// even no of dimensions
	if(delta_%2==0) {
	unsigned int n = delta_/2;
	prefactor_ = 2.pow(Constants::pi,int(n));
	for(unsigned int ix=1;ix<n;++ix) {
	prefactor_ /= double(ix);
	}
	}
	// odd number of dimensions
	else {
	unsigned int n = (delta_-1)/2;
	prefactor_ = 2.pow(Constants::pi,int(n));
	for(unsigned int ix=0;ix<n;++ix) {
	prefactor_ /= double(ix)+0.5;
	}
	}
	}

	Energy GravitonMassGenerator::mass(double & wgt, const ParticleData & ,
	const Energy low,const Energy upp, int,
	double r) const {
	Energy low2 = max(mMin_,low);
	double rlow = pow(double(low2/md_),int(delta_))/double(delta_);
	double rupp = pow(double(upp /md_),int(delta_))/double(delta_);
	double rho = rlow + (rupp-rlow)*r;
	wgt = (rupp-rlow)*prefactor_;
	return pow(double(delta_)rho,1./double(delta_))md_;
	}
	diff --git a/PDT/GenericMassGenerator.cc b/PDT/GenericMassGenerator.cc
	--- a/PDT/GenericMassGenerator.cc
	+++ b/PDT/GenericMassGenerator.cc
	@@ -1,232 +1,274 @@
	// -- C++ --
	//
	// GenericMassGenerator.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 GenericMassGenerator class.
	//

	#include "GenericMassGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/Parameter.h"
	#include "ThePEG/Interface/Switch.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "ThePEG/Utilities/Throw.h"
	#include "ThePEG/Utilities/Rebinder.h"
	+#include "GenericWidthGenerator.h"

	using namespace Herwig;
	using namespace ThePEG;

	+GenericMassGenerator::GenericMassGenerator()
	+ : maxWgt_(),BWShape_(0),nGenerate_(100),
	+ lowerMass_(),upperMass_(),
	+ mass_(),width_(),mass2_(),mWidth_(),
	+ nInitial_(1000),
	+ initialize_(false), output_(false), widthOpt_(false) {
	+}
	+
	+GenericMassGenerator::~GenericMassGenerator() {
	+}
	+
	+IBPtr GenericMassGenerator::clone() const {return new_ptr(*this);}
	+IBPtr GenericMassGenerator::fullclone() const {return new_ptr(*this);}
	+
	void GenericMassGenerator::persistentOutput(PersistentOStream & os) const {
	os << particle_
	<< ounit(lowerMass_,GeV) << ounit(upperMass_,GeV) << maxWgt_
	<< BWShape_ << nGenerate_
	<< ounit(mass_,GeV) << ounit(width_,GeV)
	<< ounit(mass2_,GeV2) << ounit(mWidth_,GeV2)
	- << nInitial_ << initialize_ << output_ << widthGen_ << widthOpt_;
	+ << nInitial_ << initialize_ << output_ << widthGen_ << widthGenB_
	+ << widthOpt_;
	}

	void GenericMassGenerator::persistentInput(PersistentIStream & is, int) {
	is >> particle_
	>> iunit(lowerMass_,GeV) >> iunit(upperMass_,GeV) >> maxWgt_
	>> BWShape_ >> nGenerate_
	>> iunit(mass_,GeV) >> iunit(width_ ,GeV)
	>> iunit(mass2_,GeV2) >> iunit(mWidth_ ,GeV2)
	- >> nInitial_ >> initialize_ >> output_ >> widthGen_ >> widthOpt_;
	+ >> nInitial_ >> initialize_ >> output_ >> widthGen_ >> widthGenB_
	+ >> widthOpt_;
	}

	ClassDescription<GenericMassGenerator> GenericMassGenerator::initGenericMassGenerator;
	// Definition of the static class description member.


	void GenericMassGenerator::Init() {

	static ClassDocumentation<GenericMassGenerator> documentation
	("The GenericMassGenerator class is the main class for"
	" mass generation in Herwig.");

	static Parameter<GenericMassGenerator,string> interfaceParticle
	("Particle",
	"The particle data object for this class",
	0, "", true, false,
	&GenericMassGenerator::setParticle,
	&GenericMassGenerator::getParticle);

	static Switch<GenericMassGenerator,bool> interfaceInitialize
	("Initialize",
	"Initialize the calculation of the maximum weight etc",
	&GenericMassGenerator::initialize_, false, false, false);
	static SwitchOption interfaceInitializeInitialization
	(interfaceInitialize,
	"Yes",
	"Do the initialization",
	true);
	static SwitchOption interfaceInitializeNoInitialization
	(interfaceInitialize,
	"No",
	"Don't do the initalization",
	false);

	static Switch<GenericMassGenerator,bool> interfaceOutput
	("Output",
	"Output the setup",
	&GenericMassGenerator::output_, false, false, false);
	static SwitchOption interfaceOutputYes
	(interfaceOutput,
	"Yes",
	"Output the data",
	true);
	static SwitchOption interfaceOutputNo
	(interfaceOutput,
	"No",
	"Don't output the data",
	false);

	static Switch<GenericMassGenerator,int> interfaceBreitWignerShape
	("BreitWignerShape",
	"Controls the shape of the mass distribution generated",
	&GenericMassGenerator::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 Parameter<GenericMassGenerator,double> interfaceMaximumWeight
	("MaximumWeight",
	"The maximum weight for the unweighting",
	&GenericMassGenerator::maxWgt_, 1.0, 0.0, 1000.0,
	false, false, true);

	static Parameter<GenericMassGenerator,int> interfaceNGenerate
	("NGenerate",
	"The number of tries to generate the mass",
	&GenericMassGenerator::nGenerate_, 100, 0, 10000,
	false, false, true);

	static Parameter<GenericMassGenerator,int> interfaceNInitial
	("NInitial",
	"Number of tries for the initialisation",
	&GenericMassGenerator::nInitial_, 1000, 0, 100000,
	false, false, true);

	static Switch<GenericMassGenerator,bool> interfaceWidthOption
	("WidthOption",
	"Which width to use",
	&GenericMassGenerator::widthOpt_, false, false, false);
	static SwitchOption interfaceWidthOptionNominalWidth
	(interfaceWidthOption,
	"NominalWidth",
	"Use the normal width from the particle data object",
	false);
	static SwitchOption interfaceWidthOptionPhysicalWidth
	(interfaceWidthOption,
	"PhysicalWidth",
	"Use the width calculated at the on-shell mass",
	true);

	}

	bool GenericMassGenerator::accept(const ParticleData & in) const {
	if(!particle_) return false;
	return in.id() == particle_->id() \|\|
	( particle_->CC() && particle_->CC()->id() == in.id() );
	}

	void GenericMassGenerator::doinit() {
	MassGenerator::doinit();
	if(particle_->massGenerator()!=this) return;
	// the width generator
	particle_->init();
	widthGen_=particle_->widthGenerator();
	- if(widthGen_){widthGen_->init();}
	+ if(widthGen_) widthGen_->init();
	+ widthGenB_ = dynamic_ptr_cast<GenericWidthGeneratorPtr>(widthGen_);
	// local storage of particle properties for speed
	mass_=particle_->mass();
	width_= widthOpt_ ? particle_->generateWidth(mass_) : particle_->width();
	mass2_=mass_*mass_;
	mWidth_=mass_*width_;
	lowerMass_ = mass_-particle_->widthLoCut();
	upperMass_ = mass_+particle_->widthUpCut();
	// print out messagw if doing the initialisation
	if(initialize_) {
	// zero the maximum weight
	maxWgt_=0.;
	// storage of variable for the loop
	double wgt=0.;
	// perform the initialisation
	// (mass() contains a ThePEG::Random call)
	for(int ix=0; ix<nInitial_; ++ix) {
	mass(wgt, *particle_, 3);
	if ( wgt > maxWgt_ )
	maxWgt_ = wgt;
	}
	}
	}

	void GenericMassGenerator::dataBaseOutput(ofstream & output, bool header) {
	if(header) output << "update Mass_Generators set parameters=\"";
	output << "newdef " << fullName() << ":BreitWignerShape " << BWShape_ << "\n";
	output << "newdef " << fullName() << ":MaximumWeight " << maxWgt_ << "\n";
	output << "newdef " << fullName() << ":NGenerate " << nGenerate_ << "\n";
	output << "newdef " << fullName() << ":WidthOption " << widthOpt_ << "\n";
	if(header) output << "\n\" where BINARY ThePEGFullName=\""
	<< fullName() << "\";" << endl;
	}

	void GenericMassGenerator::dofinish() {
	if(output_) {
	string fname = CurrentGenerator::current().filename() +
	string("-") + name() + string(".output");
	ofstream output(fname.c_str());
	dataBaseOutput(output,true);
	}
	MassGenerator::dofinish();
	}

	void GenericMassGenerator::setParticle(string p) {
	if ( (particle_ = Repository::GetPtr<tPDPtr>(p)) ) return;
	particle_ = Repository::findParticle(StringUtils::basename(p));
	if ( ! particle_ )
	Throw<InterfaceException>()
	<< "Could not set Particle interface "
	<< "for the object \"" << name()
	<< "\". Particle \"" << StringUtils::basename(p) << "\" not found."
	<< Exception::runerror;
	}

	string GenericMassGenerator::getParticle() const {
	return particle_ ? particle_->fullName() : "";
	}

	void GenericMassGenerator::rebind(const TranslationMap & trans)
	{
	particle_ = trans.translate(particle_);
	MassGenerator::rebind(trans);
	}

	IVector GenericMassGenerator::getReferences() {
	IVector ret = MassGenerator::getReferences();
	ret.push_back(particle_);
	return ret;
	}
	+
	+pair<Energy,Energy> GenericMassGenerator::width(Energy q,int shape) const {
	+ Energy gam;
	+ if(shape==3 \|\| BWShape_==1 \|\| !widthGen_ ) {
	+ gam = width_;
	+ }
	+ else if(!widthGenB_) {
	+ gam = widthGen_->width(*particle_,q);
	+ }
	+ else {
	+ return widthGenB_->width(q,*particle_);
	+ }
	+ pair<Energy,Energy> output;
	+ output.first = width_;
	+ output.second = width_;
	+ // take BR's into account
	+ double brSum=0.;
	+ for(DecaySet::const_iterator it = particle_->decayModes().begin();
	+ it!=particle_->decayModes().end();++it) {
	+ if((it).on()) brSum += (it).brat();
	+ }
	+ output.first *= brSum;
	+ return output;
	+}
	diff --git a/PDT/GenericMassGenerator.h b/PDT/GenericMassGenerator.h
	--- a/PDT/GenericMassGenerator.h
	+++ b/PDT/GenericMassGenerator.h
	@@ -1,496 +1,499 @@
	// -- C++ --
	//
	// GenericMassGenerator.h is a part of Herwig - A multi-purpose Monte Carlo event generator
	// Copyright (C) 2002-2011 The Herwig Collaboration
	//
	// Herwig is licenced under version 2 of the GPL, see COPYING for details.
	// Please respect the MCnet academic guidelines, see GUIDELINES for details.
	//
	#ifndef HERWIG_GenericMassGenerator_H
	#define HERWIG_GenericMassGenerator_H
	//
	// This is the declaration of the GenericMassGenerator class.
	//
	#include "ThePEG/PDT/MassGenerator.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/PDT/ParticleData.h"
	+#include "GenericMassGenerator.fh"
	#include "ThePEG/PDT/WidthGenerator.h"
	-#include "GenericMassGenerator.fh"
	+#include "GenericWidthGenerator.fh"
	#include "ThePEG/Repository/CurrentGenerator.h"

	namespace Herwig {
	using namespace ThePEG;

	/**
	* Declare ModelGenerator class as must be friend to set the particle
	*/
	class ModelGenerator;

	/** \ingroup PDT
	*
	* The <code>GenericMassGenerator</code> class is a simple class for the
	* generation of particle masses in Herwig. It inherits from the
	* <code>MassGenerator</code> class of ThePEG and implements a Breit-Wigner
	* using the width generator to give the running width.
	*
	* In general the width generator will be an instance of the
	* <code>GenericWidthGenerator</code> class which uses the Herwig decayers
	* based on the <code>DecayIntegrator</code> class to define the shape of the
	* running width.
	*
	* This class is designed so that the weight
	*
	* \f[\int dm^2 \frac{m\Gamma(m)}{(m^2-M^2)^2+m^2\Gamma^2(m)}\f]
	*
	* can be included in the production of the particle to take off-shell effects into
	* account. This is the default form of the weight.
	* The numerator and running of the width can
	* be changed using the BreitWignerShape interface.
	*
	* @see MassGenerator
	* @see DecayIntegrator
	* @see GenericWidthGenerator
	*
	*/
	class GenericMassGenerator: public MassGenerator {

	/**
	* ModelGenerator class as must be friend to set the particle
	*/
	friend class ModelGenerator;

	public:

	/**
	* Default constructor
	*/
	- GenericMassGenerator()
	- : maxWgt_(),BWShape_(0),nGenerate_(100),
	- lowerMass_(),upperMass_(),
	- mass_(),width_(),mass2_(),mWidth_(),
	- nInitial_(1000),
	- initialize_(false), output_(false), widthOpt_(false) {}
	+ GenericMassGenerator();
	+
	+ /**
	+ * Destructor
	+ */
	+ virtual ~GenericMassGenerator();

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

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

	/**
	* Standard Init function used to initialize the interfaces.
	*/
	static void Init();

	public:

	/**
	* Return true if this mass generator can handle the given particle type.
	* @param part The particle data pointer of the particle.
	* @return True ig this class can handle the particle and false otherwise
	*/
	bool accept(const ParticleData & part) const;

	/** @name Members to generate the mass of a particle instance */
	//@{
	/**
	* Generate a mass using the default limits.
	* @param part The particle data pointer of the particle.
	* @return The mass of the particle instance.
	*/
	Energy mass(const ParticleData & part) const {
	return mass(part,lowerMass_,upperMass_);
	}

	/**
	* Generate a mass using specified limits
	* @param part The particle data pointer of the particle.
	* @param low The lower limit on the particle's mass.
	* @param upp The upper limit on the particle's mass.
	* @return The mass of the particle instance.
	*/
	Energy mass(const ParticleData & part,
	const Energy low,const Energy upp) const {
	if(upp<low) return low;
	Energy output;
	int ntry=0; double wgt=0.;
	do {
	++ntry;
	output=mass(wgt,part,low,upp,3);
	if(wgt>maxWgt_) maxWgt_=wgt;
	}
	while(maxWgt_*(UseRandom::rnd())>wgt&&ntry<nGenerate_);
	return (ntry>=nGenerate_) ? mass_ : output;
	}

	/**
	* Return a mass with the weight using the default limits.
	* @param part The particle data pointer of the particle.
	* @param wgt The weight for this mass.
	* @param r The random number used for the weight
	* @return The mass of the particle instance.
	*/
	Energy mass(double & wgt, const ParticleData & part,
	double r=UseRandom::rnd()) const {
	return mass(wgt,part,lowerMass_,upperMass_,r);
	}

	/**
	* Return a mass with the weight using the specified limits.
	* @param part The particle data pointer of the particle.
	* @param low The lower limit on the particle's mass.
	* @param upp The upper limit on the particle's mass.
	* @param wgt The weight for this mass.
	* @param r The random number used for the weight
	* @return The mass of the particle instance.
	*/
	Energy mass(double & wgt, const ParticleData & part,
	const Energy low,const Energy upp,
	double r=UseRandom::rnd()) const {
	return mass(wgt,part,low,upp,BWShape_,r);
	}

	/**
	* Weight for the factor.
	* @param q The mass of the instance
	* @return The weight.
	*/
	virtual double weight(Energy q) const {
	return weight(q,BWShape_);
	}

	/**
	* Return the full weight
	*/
	virtual InvEnergy2 BreitWignerWeight(Energy q) {
	return BreitWignerWeight(q,BWShape_);
	}
	//@}

	/**
	* Output the initialisation info for the database
	*/
	virtual void dataBaseOutput(ofstream &,bool);

	public:

	/** @name Access to particle properties */
	//@{
	/**
	* The running width.
	* @param q The mass for the calculation of the running width
	* @return The running width.
	*/
	- Energy width(Energy q) const {
	- return (BWShape_==1\|\|!widthGen_) ?
	- width_ : widthGen_->width(*particle_,q);
	- }
	+ pair<Energy,Energy> width(Energy q,int shape) const;

	/**
	* Lower limit on the mass
	*/
	Energy lowerLimit() const {return lowerMass_;}

	/**
	* Upper limit on the mass
	*/
	Energy upperLimit() const {return upperMass_;}

	/**
	* Default mass
	*/
	Energy nominalMass() const {return mass_;}

	/**
	* Default Width
	*/
	Energy nominalWidth() const {return width_;}

	protected:

	/**
	* Return a mass with the weight using the specified limits.
	* @param low The lower limit on the particle's mass.
	* @param upp The upper limit on the particle's mass.
	* @param wgt The weight for this mass.
	* @param shape The type of shape to use
	* @param r The random number used for the weight
	* @return The mass of the particle instance.
	*/
	virtual Energy mass(double & wgt, const ParticleData & ,
	const Energy low,const Energy upp, int shape,
	double r=UseRandom::rnd()) const {
	// calculate the mass square using fixed width BW
	Energy lo=max(low,lowerMass_),up=min(upp,upperMass_);
	double rhomin=atan2((lo*lo-mass2_),mWidth_);
	double rhomax=atan2((up*up-mass2_),mWidth_)-rhomin;
	double rho=rhomin+rhomax*r;
	Energy2 q2 = mass2_+mWidth_*tan(rho);
	Energy q = sqrt(q2);
	wgt = rhomax*weight(q,shape);
	// return the mass
	return q;
	}

	/**
	* Return a mass with the weight using the default limits.
	* @param part The particle data pointer of the particle.
	* @param wgt The weight for this mass.
	* @param shape The type of shape to use
	* @param r The random number used for the weight
	* @return The mass of the particle instance.
	*/
	Energy mass(double & wgt, const ParticleData & part, int shape,
	double r=UseRandom::rnd()) const {
	return mass(wgt,part,lowerMass_,upperMass_,shape,r);
	}

	/**
	* Weight for the factor.
	* @param q The mass of the instance
	* @param shape The type of shape to use as for the BreitWignerShape interface
	* @return The weight.
	*/
	inline virtual double weight(Energy q, int shape) const {
	Energy2 q2 = q*q;
	Energy4 sq=sqr(q2-mass2_);
	- Energy gam=width(q);
	+ pair<Energy,Energy> gam=width(q,shape);
	// finish the calculation of the width
	Energy2 num;
	- if(shape==2) num = mass_*gam;
	- else if(shape==3) num = mass_*width_;
	- else num = q*gam;
	- Energy4 den = (shape==2) ? sq+mass2_gamgam : sq+q2gamgam;
	+ if(shape==2) num = mass_*gam.first;
	+ else if(shape==3) num = mass_*gam.first;
	+ else num = q *gam.first;
	+ Energy4 den = (shape==2) ? sq+mass2_gam.secondgam.second : sq+q2gam.secondgam.second;
	return num/den(sq+mWidth_mWidth_)/Constants::pi/mWidth_;
	}

	/**
	* Return the full weight
	*/
	virtual InvEnergy2 BreitWignerWeight(Energy q, int shape) const {
	Energy2 q2 = q*q;
	Energy4 sq=sqr(q2-mass2_);
	- Energy gam=width(q);
	+ pair<Energy,Energy> gam=width(q,shape);
	// finish the calculation of the width
	Energy2 num;
	- if(shape==2) num = mass_*gam;
	- else if(shape==3) num = mass_*width_;
	- else num = q*gam;
	- Energy4 den = (shape==2) ? sq+mass2_gamgam : sq+q2gamgam;
	+ if(shape==2) num = mass_*gam.first;
	+ else if(shape==3) num = mass_*gam.first;
	+ else num = q *gam.first;
	+ Energy4 den = (shape==2) ? sq+mass2_gam.secondgam.second : sq+q2gam.secondgam.second;
	return num/den/Constants::pi;
	}

	/**
	* Accesss to the particle
	*/
	tcPDPtr particle() const {return particle_;}

	/**
	* Set the particle
	*/
	void particle(tPDPtr in) {particle_ = in;}

	protected:

	/** @name Clone Methods. */
	//@{
	/**
	* Make a simple clone of this object.
	* @return a pointer to the new object.
	*/
	- virtual IBPtr clone() const {return new_ptr(*this);}
	+ 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 {return new_ptr(*this);}
	+ virtual IBPtr fullclone() const;
	//@}

	protected:

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

	/**
	* Finalize this object. Called in the run phase just after a
	* run has ended. Used eg. to write out statistics.
	*/
	virtual void dofinish();

	/**
	* Rebind pointer to other Interfaced objects. Called in the setup phase
	* after all objects used in an EventGenerator has been cloned so that
	* the pointers will refer to the cloned objects afterwards.
	* @param trans a TranslationMap relating the original objects to
	* their respective clones.
	* @throws RebindException if no cloned object was found for a given
	* pointer.
	*/
	virtual void rebind(const TranslationMap & trans)
	;

	/**
	* Return a vector of all pointers to Interfaced objects used in this
	* object.
	* @return a vector of pointers.
	*/
	virtual IVector getReferences();
	//@}

	private:

	/**
	* Describe a concrete class with persistent data.
	*/
	static ClassDescription<GenericMassGenerator> initGenericMassGenerator;

	/**
	* Private and non-existent assignment operator.
	*/
	GenericMassGenerator & operator=(const GenericMassGenerator &);

	private:

	/**
	* Helper function for the interface
	*/
	void setParticle(string);

	/**
	* Helper function for the interface
	*/
	string getParticle() const;

	private:

	/**
	* The maximum weight for unweighting when generating the mass.
	*/
	mutable double maxWgt_;

	/**
	* parameter controlling the shape of the Breit-Wigner
	*/
	int BWShape_;

	/**
	* Number of attempts to generate the mass.
	*/
	int nGenerate_;

	private:

	/**
	* Pointer to the particle
	*/
	tPDPtr particle_;

	/**
	* Lower limit on the particle's mass
	*/
	Energy lowerMass_;

	/**
	* Upper limit on the particle's mass
	*/
	Energy upperMass_;

	/**
	* Mass of the particle
	*/
	Energy mass_;

	/**
	* Width of the particle
	*/
	Energy width_;

	/**
	* Mass of the particle squared.
	*/
	Energy2 mass2_;

	/**
	* Mass of the particle times the width.
	*/
	Energy2 mWidth_;

	/**
	* Number of weights to generate when initializing
	*/
	int nInitial_;

	/**
	* Whether or not to initialize the GenericMassGenerator
	*/
	bool initialize_;

	/**
	* Whether or not to output the data to a file
	*/
	bool output_;

	/**
	* Pointer to the width generator
	*/
	WidthGeneratorPtr widthGen_;

	/**
	+ * Pointer to the width generator
	+ */
	+ GenericWidthGeneratorPtr widthGenB_;
	+
	+ /**
	* Option for the treatment of the width
	*/
	bool widthOpt_;

	};

	}


	#include "ThePEG/Utilities/ClassTraits.h"

	namespace ThePEG {

	/** @cond TRAITSPECIALIZATIONS */

	/**
	* The following template specialization informs ThePEG about the
	* base class of GenericMassGenerator.
	*/
	template <>
	struct BaseClassTrait<Herwig::GenericMassGenerator,1> {
	/** Typedef of the base class of GenericMassGenerator. */
	typedef MassGenerator NthBase;
	};

	/**
	* The following template specialization informs ThePEG about the
	* name of this class and the shared object where it is defined.
	*/
	template <>
	struct ClassTraits<Herwig::GenericMassGenerator>
	: public ClassTraitsBase<Herwig::GenericMassGenerator> {
	/** Return the class name. */
	static string className() { return "Herwig::GenericMassGenerator"; }
	};

	/** @endcond */

	}

	#endif /* HERWIG_GenericMassGenerator_H */
	diff --git a/PDT/ScalarMassGenerator.cc b/PDT/ScalarMassGenerator.cc
	--- a/PDT/ScalarMassGenerator.cc
	+++ b/PDT/ScalarMassGenerator.cc
	@@ -1,163 +1,164 @@
	// -- C++ --
	//
	// ScalarMassGenerator.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 ScalarMassGenerator class.
	//

	#include "ScalarMassGenerator.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/Interface/ParVector.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	+#include "GenericWidthGenerator.h"

	using namespace Herwig;

	void ScalarMassGenerator::persistentOutput(PersistentOStream & os) const {
	os << ounit(_coupling,GeV) << ounit(_mass1,GeV) << ounit(_mass2,GeV)
	<< ounit(_m2plus,GeV2) << ounit(_m2minus,GeV2);
	}

	void ScalarMassGenerator::persistentInput(PersistentIStream & is, int) {
	is >> iunit(_coupling,GeV) >> iunit(_mass1,GeV) >> iunit(_mass2,GeV)
	>> iunit(_m2plus,GeV2) >> iunit(_m2minus,GeV2);
	}

	ClassDescription<ScalarMassGenerator> ScalarMassGenerator::initScalarMassGenerator;
	// Definition of the static class description member.

	void ScalarMassGenerator::Init() {

	static ClassDocumentation<ScalarMassGenerator> documentation
	("The ScalarMassGenerator class inherits from the "
	"GenericMassGenerator class and includes finite width effects "
	"for the scalar f_0 and a_0 mesons.",
	"Finite width effects for the scalar $f_0$ and $a_0$ mesons follow \\cite{Flatte:1976xu}.",
	"%\\cite{Flatte:1976xu}\n"
	"\\bibitem{Flatte:1976xu}\n"
	" S.~M.~Flatte,\n"
	" ``Coupled - Channel Analysis Of The Pi Eta And K Anti-K Systems Near K Anti-K\n"
	" Threshold,''\n"
	" Phys.\\ Lett.\\ B {\\bf 63}, 224 (1976).\n"
	" %%CITATION = PHLTA,B63,224;%%\n"
	);

	static ParVector<ScalarMassGenerator,Energy> interfacecoupling
	("Coupling",
	"The coupling",
	&ScalarMassGenerator::_coupling,
	GeV, 0, ZERO, ZERO, Constants::MaxEnergy, false, false, true);

	static ParVector<ScalarMassGenerator,Energy> interfacemass1
	("Mass1",
	"The mass for first particle",
	&ScalarMassGenerator::_mass1,
	GeV, 0, ZERO, ZERO, Constants::MaxEnergy, false, false, true);

	static ParVector<ScalarMassGenerator,Energy> interfacemass2
	("Mass2",
	"The mass for second particle",
	&ScalarMassGenerator::_mass2,
	GeV, 0, ZERO, ZERO, Constants::MaxEnergy, false, false, true);

	}

	void ScalarMassGenerator::dataBaseOutput(ofstream & output,bool header) {
	if(header) output << "update Mass_Generators set parameters=\"";
	GenericMassGenerator::dataBaseOutput(output,false);
	for(unsigned int ix=0;ix<_coupling.size();++ix)
	output << "insert " << fullName() << ":Coupling "
	<< ix << " " << _coupling[ix]/GeV << "\n";
	for(unsigned int ix=0;ix<_mass1.size();++ix)
	output << "insert " << fullName() << ":Mass1 "
	<< ix << " " << _mass1[ix]/GeV << "\n";
	for(unsigned int ix=0;ix<_mass2.size();++ix)
	output << "insert " << fullName()
	<< ":Mass2 " << ix << " " << _mass2[ix]/GeV << "\n";
	if(header) output << "\n\" where BINARY ThePEGName=\""
	<< fullName() << "\";" << endl;
	}

	void ScalarMassGenerator::doinit() {
	if(_coupling.size()!=_mass1.size()\|\|
	_coupling.size()!=_mass2.size())
	throw InitException() << "Parameter vectors have inconsistent sizes in "
	<< "ScalarMassGenerator::doinit()"
	<< Exception::runerror;
	// initialise the local variables
	for(unsigned int ix=0;ix<_mass1.size();++ix) {
	_m2plus .push_back(sqr(_mass1[ix]+_mass2[ix]));
	_m2minus.push_back(sqr(_mass1[ix]-_mass2[ix]));
	}
	// rest of the initialisation is handled in the base class
	GenericMassGenerator::doinit();
	}

	double ScalarMassGenerator::weight(Energy q,int shape) const {
	useMe();
	Energy2 q2 = sqr(q);
	Energy2 mass2 = sqr(nominalMass());
	Energy2 mwidth= nominalMass()*nominalWidth();
	Energy2 gamma[2]={0.*MeV2,ZERO};
	if(shape==1) {
	for(unsigned int ix=0;ix<_coupling.size();++ix) {
	double lambda = (mass2-_m2plus[ix])*(mass2-_m2minus[ix])/sqr(mass2);
	if(lambda>=0.) gamma[0]+=sqr(_coupling[ix])sqrt( lambda)q/nominalMass();
	else gamma[1]+=sqr(_coupling[ix])sqrt(-lambda)q/nominalMass();
	}
	}
	else {
	for(unsigned int ix=0;ix<_coupling.size();++ix) {
	double lambda = (q2-_m2plus[ix])*(q2-_m2minus[ix])/sqr(q2);
	if(lambda>=0.) gamma[0]+=sqr(_coupling[ix])*sqrt( lambda);
	else gamma[1]+=sqr(_coupling[ix])*sqrt(-lambda);
	}
	}
	Energy2 numer(ZERO);
	if(shape==2) numer = nominalMass()*gamma[0]/q;
	else if(shape==3) numer = nominalMass()*nominalWidth();
	else numer = gamma[0];
	complex<Energy2> denom = (shape==2) ?
	mass2 - q2 +mass2/q2*complex<Energy2>(gamma[1],-gamma[0]) :
	mass2 - q2 +complex<Energy2>(gamma[1],-gamma[0]);
	// complex denominantor
	Energy4 den = real(denom*conj(denom));
	return numer/den*(sqr(mass2-q2)+sqr(mwidth))/Constants::pi/mwidth;
	}

	InvEnergy2 ScalarMassGenerator::BreitWignerWeight(Energy q, int shape) const {
	useMe();
	Energy2 q2 = sqr(q);
	Energy2 mass2 = sqr(nominalMass());
	Energy2 gamma[2]={0.*MeV2,ZERO};
	if(shape==1) {
	for(unsigned int ix=0;ix<_coupling.size();++ix) {
	double lambda = (mass2-_m2plus[ix])*(mass2-_m2minus[ix])/sqr(mass2);
	if(lambda>=0.) gamma[0]+=sqr(_coupling[ix])sqrt( lambda)q/nominalMass();
	else gamma[1]+=sqr(_coupling[ix])sqrt(-lambda)q/nominalMass();
	}
	}
	else {
	for(unsigned int ix=0;ix<_coupling.size();++ix) {
	double lambda = (q2-_m2plus[ix])*(q2-_m2minus[ix])/sqr(q2);
	if(lambda>=0.) gamma[0]+=sqr(_coupling[ix])*sqrt( lambda);
	else gamma[1]+=sqr(_coupling[ix])*sqrt(-lambda);
	}
	}
	Energy2 numer(ZERO);
	if(shape==2) numer = nominalMass()*gamma[0]/q;
	else if(shape==3) numer = nominalMass()*nominalWidth();
	else numer = gamma[0];
	complex<Energy2> denom = (shape==2) ?
	mass2 - q2 +mass2/q2*complex<Energy2>(gamma[1],-gamma[0]) :
	mass2 - q2 +complex<Energy2>(gamma[1],-gamma[0]);
	// complex denominantor
	Energy4 den = real(denom*conj(denom));
	return numer/den/Constants::pi;
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 7:51 PM (1 d, 9 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3792764
Default Alt Text: (33 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions