No OneTemporary
Actions

Size

28 KB

Subscribers

None

View Options

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

	#include "PhaseSpaceMode.h"
	#include "ThePEG/Interface/ClassDocumentation.h"
	#include "ThePEG/EventRecord/Particle.h"
	#include "ThePEG/Repository/UseRandom.h"
	#include "ThePEG/Repository/EventGenerator.h"
	#include "ThePEG/Utilities/DescribeClass.h"
	#include "ThePEG/Persistency/PersistentOStream.h"
	#include "ThePEG/Persistency/PersistentIStream.h"
	#include "Herwig/Utilities/Kinematics.h"

	using namespace Herwig;
	using namespace ThePEG::Helicity;

	void PhaseSpaceMode::persistentOutput(PersistentOStream & os) const {
	os << incoming_ << channels_ << maxWeight_ << outgoing_ << outgoingCC_
	<< partial_ << widthGen_ << massGen_ << testOnShell_
	<< ounit(eMax_,GeV) << nRand_;
	}

	void PhaseSpaceMode::persistentInput(PersistentIStream & is, int) {
	is >> incoming_ >> channels_ >> maxWeight_ >> outgoing_ >> outgoingCC_
	>> partial_ >> widthGen_ >> massGen_ >> testOnShell_
	>> iunit(eMax_,GeV) >> nRand_;
	}

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

	void PhaseSpaceMode::Init() {

	static ClassDocumentation<PhaseSpaceMode> documentation
	("There is no documentation for the PhaseSpaceMode class");

	}

	ParticleVector
	PhaseSpaceMode::generateDecay(const Particle & inpart,
	tcDecayIntegrator2Ptr decayer,
	bool intermediates,bool cc) {
	decayer->ME(DecayMEPtr());
	eps_=decayer->eps_;
	// compute the prefactor
	Energy prewid = (widthGen_&&partial_>=0) ?
	widthGen_->partialWidth(partial_,inpart.mass()) :
	incoming_.first->width();
	InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
	// boosts to/from rest
	Boost bv =-inpart.momentum().boostVector();
	double gammarest = inpart.momentum().e()/inpart.momentum().mass();
	LorentzRotation boostToRest( bv,gammarest);
	LorentzRotation boostFromRest(-bv,gammarest);
	// construct a new particle which is at rest
	Particle inrest(inpart);
	inrest.transform(boostToRest);
	double wgt(0.);
	vector<Lorentz5Momentum> momenta(outgoing_.size());
	int ichan;
	unsigned int ncount(0);
	try {
	do {
	// phase-space pieces of the weight
	fillStack();
	- wgt = pre*weight(ichan,inrest,momenta);
	+ wgt = pre*weight(ichan,inrest.momentum(),momenta);
	// matrix element piece
	wgt *= decayer->me2(-1,inrest,!cc ? outgoing_ : outgoingCC_,
	momenta,
	ncount ==0 ? DecayIntegrator2::Initialize : DecayIntegrator2::Calculate);
	++ncount;
	if(wgt>maxWeight_) {
	CurrentGenerator::log() << "Resetting max weight for decay "
	<< inrest.PDGName() << " -> ";
	for(tcPDPtr part : outgoing_)
	CurrentGenerator::log() << " " << part->PDGName();
	CurrentGenerator::log() << " " << maxWeight_ << " " << wgt
	<< " " << inrest.mass()/MeV << "\n";
	maxWeight_=wgt;
	}
	}
	while(maxWeight_*UseRandom::rnd()>wgt&&ncount<decayer->nTry_);
	}
	catch (Veto) {
	// restore the incoming particle to its original state
	inrest.transform(boostFromRest);
	throw Veto();
	}
	if(ncount>=decayer->nTry_) {
	CurrentGenerator::log() << "The decay " << inrest.PDGName() << " -> ";
	for(tcPDPtr part : outgoing_)
	CurrentGenerator::log() << " " << part->PDGName();
	CurrentGenerator::log() << " " << maxWeight_ << " " << decayer->nTry_
	<< " is too inefficient for the particle "
	<< inpart << "vetoing the decay \n";
	momenta.clear();
	throw Veto();
	}
	// construct the particles
	ParticleVector output(outgoing_.size());
	for(unsigned int ix=0;ix<outgoing_.size();++ix)
	output[ix] = (!cc ? outgoing_[ix] : outgoingCC_[ix] )->produceParticle(momenta[ix]);
	// set up the spin correlations
	decayer->constructSpinInfo(inrest,output);
	const_ptr_cast<tPPtr>(&inpart)->spinInfo(inrest.spinInfo());
	constructVertex(inpart,output,decayer);
	// return if intermediate particles not required
	if(channels_.empty()\|\|!intermediates) {
	for(tPPtr part : output) part->transform(boostFromRest);
	}
	// find the intermediate particles
	else {
	// select the channel
	vector<double> mewgts(channels_.size(),0.0);
	double total=0.;
	for(unsigned int ix=0,N=channels_.size();ix<N;++ix) {
	mewgts[ix]=decayer->me2(ix,inrest,!cc ? outgoing_ : outgoingCC_,
	momenta,DecayIntegrator2::Calculate);
	total+=mewgts[ix];
	}
	// randomly pick a channel
	total *= UseRandom::rnd();
	int iChannel = -1;
	do {
	++iChannel;
	total-=mewgts[iChannel];
	}
	while(iChannel<int(channels_.size()) && total>0.);
	// apply boost
	for(tPPtr part : output) part->transform(boostFromRest);
	channels_[iChannel].generateIntermediates(cc,inpart,output);
	}
	decayer->ME(DecayMEPtr());
	// return particles;
	return output;
	}

	// flat phase space generation and weight
	-Energy PhaseSpaceMode::flatPhaseSpace(const Particle & inpart,
	+Energy PhaseSpaceMode::flatPhaseSpace(const Lorentz5Momentum & in,
	vector<Lorentz5Momentum> & momenta,
	bool onShell) const {
	double wgt(1.);
	// masses of the particles
	- vector<Energy> mass = externalMasses(inpart.mass(),wgt,onShell);
	+ vector<Energy> mass = externalMasses(in.mass(),wgt,onShell);
	// two body decay
	double ctheta = 2.*rStack_.top() - 1.;
	rStack_.pop();
	double phi = Constants::twopi*rStack_.top();
	rStack_.pop();
	- if(! Kinematics::twoBodyDecay(inpart.momentum(), mass[0], mass[1],
	+ if(! Kinematics::twoBodyDecay(in, mass[0], mass[1],
	ctheta, phi, momenta[0],momenta[1]))
	throw Exception() << "Incoming mass - Outgoing mass negative in "
	<< "PhaseSpaceMode::flatPhaseSpace()"
	<< Exception::eventerror;
	- wgt *= Kinematics::pstarTwoBodyDecay(inpart.mass(),mass[0],mass[1])
	- /8./Constants::pi/inpart.mass();
	- return wgt*inpart.mass();
	+ wgt *= Kinematics::pstarTwoBodyDecay(in.mass(),mass[0],mass[1])
	+ /8./Constants::pi/in.mass();
	+ return wgt*in.mass();
	}

	// generate the masses of the external particles
	vector<Energy> PhaseSpaceMode::externalMasses(Energy inmass,double & wgt,
	bool onShell) const {
	vector<Energy> mass(outgoing_.size());
	vector<int> notdone;
	Energy mlow(ZERO);
	// set masses of stable particles and limits
	for(unsigned int ix=0;ix<outgoing_.size();++ix) {
	// get the mass of the particle if can't use weight
	if(onShell)
	mass[ix] = outgoing_[ix]->mass();
	else if(!massGen_[ix] \|\| outgoing_[ix]->stable()) {
	mass[ix] = outgoing_[ix]->generateMass();
	mlow += mass[ix];
	}
	else {
	mass[ix] = ZERO;
	notdone.push_back(ix);
	mlow+=max(outgoing_[ix]->mass()-outgoing_[ix]->widthLoCut(),ZERO);
	}
	}
	if(mlow>inmass) throw Veto();
	// now we need to generate the masses for the particles we haven't
	for( ;!notdone.empty();) {
	unsigned int iloc=long(UseRandom::rnd()*(notdone.size()-1));
	Energy low = max(outgoing_[notdone[iloc]]->mass()-outgoing_[notdone[iloc]]->widthLoCut(),ZERO);
	mlow-=low;
	double wgttemp;
	mass[notdone[iloc]]= massGen_[notdone[iloc]]->mass(wgttemp,*outgoing_[notdone[iloc]],low,inmass-mlow,rStack_.top());
	rStack_.pop();
	assert(mass[notdone[iloc]]>=low&&mass[notdone[iloc]]<=inmass-mlow);
	wgt *= wgttemp;
	mlow += mass[notdone[iloc]];
	notdone.erase(notdone.begin()+iloc);
	}
	return mass;
	}

	// construct the vertex for spin corrections
	void PhaseSpaceMode::constructVertex(const Particle & inpart,
	const ParticleVector & decay,
	tcDecayIntegrator2Ptr decayer) const {
	// construct the decay vertex
	VertexPtr vertex(new_ptr(DecayVertex()));
	DVertexPtr Dvertex(dynamic_ptr_cast<DVertexPtr>(vertex));
	// set the incoming particle for the decay vertex
	(inpart.spinInfo())->decayVertex(vertex);
	for(unsigned int ix=0;ix<decay.size();++ix) {
	decay[ix]->spinInfo()->productionVertex(vertex);
	}
	// set the matrix element
	Dvertex->ME(decayer->ME());
	decayer->ME(DecayMEPtr());
	}

	// initialise the phase space
	Energy PhaseSpaceMode::initializePhaseSpace(bool init, tcDecayIntegrator2Ptr decayer,
	bool onShell) {
	decayer->ME(DecayMEPtr());
	eps_=decayer->eps_;
	Energy output(ZERO);
	// ensure that the weights add up to one
	if(!channels_.empty()) {
	double temp=0.;
	for(const PhaseSpaceChannel & channel : channels_) temp+= channel.weight();
	for(PhaseSpaceChannel & channel : channels_) channel.weight(channel.weight()/temp);
	}
	if(!init) return ZERO;
	ThePEG::PPtr inpart=incoming_.first->produceParticle();
	// now if using flat phase space
	maxWeight_=0.;
	if(channels_.empty()) {
	vector<Lorentz5Momentum> momenta(outgoing_.size());
	double wsum=0.,wsqsum=0.;
	Energy m0,mmin(ZERO);
	for(tcPDPtr part : outgoing_) mmin += part->massMin();
	for(unsigned int ix=0;ix<decayer->nPoint_;++ix) {
	// set the mass of the decaying particle
	m0 = !onShell ? inpart->dataPtr()->generateMass() : inpart->dataPtr()->mass();
	double wgt=0.;
	if(m0<=mmin) continue;
	inpart->set5Momentum(Lorentz5Momentum(m0));
	// compute the prefactor
	Energy prewid = (widthGen_&&partial_>=0) ?
	widthGen_->partialWidth(partial_,inpart->mass()) :
	incoming_.first->width();
	InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
	// generate the weight for this point
	try {
	int dummy;
	// phase-space piece
	fillStack();
	- wgt = preweight(dummy,inpart,momenta,onShell);
	+ wgt = pre*weight(dummy,inpart->momentum(),momenta,onShell);
	// matrix element piece
	wgt = decayer->me2(-1,inpart,outgoing_,momenta,DecayIntegrator2::Initialize);
	}
	catch (Veto) {
	wgt=0.;
	}
	if(wgt>maxWeight_) maxWeight_ = wgt;
	wsum += wgt;
	wsqsum += sqr(wgt);
	}
	wsum /= decayer->nPoint_;
	wsqsum=wsqsum/decayer->nPoint_-sqr(wsum);
	if(wsqsum<0.) wsqsum=0.;
	wsqsum=sqrt(wsqsum/decayer->nPoint_);
	Energy fact = (widthGen_&&partial_>=0) ?
	widthGen_->partialWidth(partial_,inpart->nominalMass()) :
	inpart->dataPtr()->width();
	if(fact==ZERO) fact=MeV;
	// factor for the weight with spin correlations
	maxWeight_ *= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6;
	if ( Debug::level > 1 ) {
	// ouptut the information on the initialisation
	CurrentGenerator::log() << "Initialized the phase space for the decay "
	<< incoming_.first->PDGName() << " -> ";
	for(tPDPtr part : outgoing_)
	CurrentGenerator::log() << part->PDGName() << " ";
	CurrentGenerator::log() << "\n";
	if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << wsum
	<< " +/- " << wsqsum << "\n";
	CurrentGenerator::log() << "The partial width is " << wsum*fact/MeV
	<< " +/- " << wsqsum*fact/MeV << " MeV\n";
	CurrentGenerator::log() << "The partial width is "
	<< wsum*fact/6.58212E-22/MeV
	<< " +/- " << wsqsum*fact/6.58212E-22/MeV<< " s-1\n";
	CurrentGenerator::log() << "The maximum weight is "
	<< maxWeight_ << endl;
	}
	output=wsum*fact;
	}
	else {
	vector<Lorentz5Momentum> momenta(outgoing_.size());
	double totsum(0.),totsq(0.);
	for(unsigned int iy=0;iy<decayer->nIter_;++iy) {
	// zero the maximum weight
	maxWeight_=0.;
	vector<double> wsum(channels_.size(),0.),wsqsum(channels_.size(),0.);
	vector<int> nchan(channels_.size(),0);
	totsum = 0.;
	totsq = 0.;
	Energy m0,mmin(ZERO);
	for(tcPDPtr part : outgoing_) mmin += part->massMin();
	for(unsigned int ix=0;ix<decayer->nPoint_;++ix) {
	m0 = !onShell ? incoming_.first->generateMass() : incoming_.first->mass();
	double wgt=0.;
	int ichan(-1);
	if(m0>mmin) {
	inpart->set5Momentum(Lorentz5Momentum(m0));
	// compute the prefactor
	Energy prewid= (widthGen_&&partial_>=0) ?
	widthGen_->partialWidth(partial_,inpart->mass()) :
	inpart->dataPtr()->width();
	InvEnergy pre = prewid>ZERO ? 1./prewid : 1./MeV;
	// generate the weight for this point
	try {
	fillStack();
	- wgt = preweight(ichan,inpart,momenta,onShell);
	+ wgt = pre*weight(ichan,inpart->momentum(),momenta,onShell);
	// matrix element piece
	wgt = decayer->me2(-1,inpart,outgoing_,momenta,DecayIntegrator2::Initialize);
	}
	catch (Veto) {
	wgt=0.;
	}
	}
	if(wgt>maxWeight_) maxWeight_=wgt;
	if(ichan>=0) {
	wsum[ichan] += wgt;
	wsqsum[ichan] += sqr(wgt);
	++nchan[ichan];
	}
	totsum+=wgt;
	totsq+=wgt*wgt;
	}
	totsum=totsum/decayer->nPoint_;
	totsq=totsq/decayer->nPoint_-sqr(totsum);
	if(totsq<0.) totsq=0.;
	totsq=sqrt(totsq/decayer->nPoint_);
	if ( Debug::level > 1 )
	CurrentGenerator::log() << "The branching ratio is " << iy << " "
	<< totsum << " +/- " << totsq
	<< maxWeight_ << "\n";
	// compute the individual terms
	double total(0.);
	for(unsigned int ix=0;ix<channels_.size();++ix) {
	if(nchan[ix]!=0) {
	wsum[ix]=wsum[ix]/nchan[ix];
	wsqsum[ix]=wsqsum[ix]/nchan[ix];
	if(wsqsum[ix]<0.) wsqsum[ix]=0.;
	wsqsum[ix]=sqrt(wsqsum[ix]/nchan[ix]);
	}
	else {
	wsum[ix]=0;
	wsqsum[ix]=0;
	}
	total+=sqrt(wsqsum[ix])*channels_[ix].weight();
	}
	if(total>0.) {
	for(unsigned int ix=0;ix<channels_.size();++ix) {
	channels_[ix].weight(sqrt(wsqsum[ix])*channels_[ix].weight()/total);
	}
	}
	}
	// factor for the weight with spin correlations
	maxWeight_*= inpart->dataPtr()->iSpin()==1 ? 1.1 : 1.6;
	// output the information on the initialisation
	Energy fact = (widthGen_&&partial_>=0) ?
	widthGen_->partialWidth(partial_,inpart->nominalMass()) :
	inpart->dataPtr()->width();
	output=totsum*fact;
	if(fact==ZERO) fact=MeV;
	if ( Debug::level > 1 ) {
	CurrentGenerator::log() << "Initialized the phase space for the decay "
	<< incoming_.first->PDGName() << " -> ";
	for(tcPDPtr part : outgoing_)
	CurrentGenerator::log() << part->PDGName() << " ";
	CurrentGenerator::log() << "\n";
	if(fact!=MeV) CurrentGenerator::log() << "The branching ratio is " << totsum
	<< " +/- " << totsq << "\n";
	CurrentGenerator::log() << "The partial width is " << totsum*fact/MeV
	<< " +/- " << totsq*fact/MeV << " MeV\n";
	CurrentGenerator::log() << "The partial width is "
	<< totsum*fact/6.58212E-22/MeV
	<< " +/- " << totsq*fact/6.58212E-22/MeV
	<< " s-1\n";
	CurrentGenerator::log() << "The maximum weight is "
	<< maxWeight_ << "\n";
	CurrentGenerator::log() << "The weights for the different phase"
	<< " space channels are \n";
	for(unsigned int ix=0,N=channels_.size();ix<N;++ix) {
	CurrentGenerator::log() << "Channel " << ix
	<< " had weight " << channels_[ix].weight()
	<< "\n";
	}
	}
	CurrentGenerator::log() << flush;
	}
	decayer->ME(DecayMEPtr());
	return output;
	}

	// generate a phase-space point using multichannel phase space
	-Energy PhaseSpaceMode::channelPhaseSpace(int & ichan, const Particle & in,
	+Energy PhaseSpaceMode::channelPhaseSpace(int & ichan, const Lorentz5Momentum & in,
	vector<Lorentz5Momentum> & momenta,
	bool onShell) const {
	double wgt(rStack_.top());
	rStack_.pop();
	// select a channel
	ichan=-1;
	do {
	++ichan;
	wgt-=channels_[ichan].weight();
	}
	while(ichan<int(channels_.size())&&wgt>0.);
	// generate the momenta
	if(ichan==int(channels_.size())) {
	throw DecayIntegratorError() << "PhaseSpaceMode::channelPhaseSpace()"
	<< " failed to select a channel"
	<< Exception::abortnow;
	}
	// generate the masses of the external particles
	double masswgt(1.);
	vector<Energy> mass(externalMasses(in.mass(),masswgt,onShell));
	- momenta=channels_[ichan].generateMomenta(in.momentum(),mass);
	+ momenta=channels_[ichan].generateMomenta(in,mass);
	// compute the denominator of the weight
	wgt=0.;
	for(const PhaseSpaceChannel & channel : channels_)
	wgt += channel.generateWeight(momenta);
	// return the weight
	return wgt!=0. ? in.mass()*masswgt/wgt : ZERO;
	}

	void PhaseSpaceMode::init() {
	// get mass and width generators
	outgoingCC_.clear();
	for(tPDPtr part : outgoing_) {
	outgoingCC_.push_back(part->CC() ? part->CC() : part);
	}
	massGen_.resize(outgoing_.size());
	widthGen_ = dynamic_ptr_cast<cGenericWidthGeneratorPtr>(incoming_.first->widthGenerator());
	for(unsigned int ix=0;ix<outgoing_.size();++ix) {
	assert(outgoing_[ix]);
	massGen_[ix]= dynamic_ptr_cast<cGenericMassGeneratorPtr>(outgoing_[ix]->massGenerator());
	}
	// get max energy for decays
	if(!incoming_.second)
	eMax_ = testOnShell_ ? incoming_.first->mass() : incoming_.first->massMax();
	// work out how many random numbers we need
	nRand_ = 3*outgoing_.size()-4;
	for(unsigned int ix=0;ix<outgoing_.size();++ix) {
	if(massGen_[ix]) ++nRand_;
	}
	if(!channels_.empty()) ++nRand_;
	}

	void PhaseSpaceMode::initrun() {
	// update the mass and width generators
	if(incoming_.first->widthGenerator()!=widthGen_)
	widthGen_ = dynamic_ptr_cast<cGenericWidthGeneratorPtr>(incoming_.first->widthGenerator());
	for(unsigned int ix=0;ix<outgoing_.size();++ix) {
	if(massGen_[ix]!=outgoing_[ix]->massGenerator())
	massGen_[ix] = dynamic_ptr_cast<cGenericMassGeneratorPtr>(outgoing_[ix]->massGenerator());
	}
	for(PhaseSpaceChannel & channel : channels_) channel.initrun(this);
	if(widthGen_) const_ptr_cast<GenericWidthGeneratorPtr>(widthGen_)->initrun();
	tcGenericWidthGeneratorPtr wtemp;
	for(unsigned int ix=0;ix<outgoing_.size();++ix) {
	wtemp=
	dynamic_ptr_cast<tcGenericWidthGeneratorPtr>(outgoing_[ix]->widthGenerator());
	if(wtemp) const_ptr_cast<tGenericWidthGeneratorPtr>(wtemp)->initrun();
	}
	}
	diff --git a/Decay/PhaseSpaceMode.h b/Decay/PhaseSpaceMode.h
	--- a/Decay/PhaseSpaceMode.h
	+++ b/Decay/PhaseSpaceMode.h
	@@ -1,372 +1,374 @@
	// -- C++ --
	#ifndef Herwig_PhaseSpaceMode_H
	#define Herwig_PhaseSpaceMode_H
	//
	// This is the declaration of the PhaseSpaceMode class.
	//

	#include "ThePEG/Config/ThePEG.h"
	#include "PhaseSpaceMode.fh"
	#include "PhaseSpaceChannel.h"
	#include "Herwig/PDT/GenericWidthGenerator.h"
	#include "Herwig/PDT/GenericMassGenerator.h"

	namespace Herwig {
	using namespace ThePEG;

	/** \ingroup Decay
	*
	* The <code>PhaseSpaceMode</code> class is designed to store a group
	* of phase-space channels for use by the DecayIntegrator class to
	* generate the phase-space for a given decay mode.
	*
	* Additional phase-space channels can be added using the addChannel member.
	*
	* In practice the modes are usually constructed together with the a number of
	* <code>PhaseSpaceChannel</code> objects. In classes inheriting from the
	* DecayIntegrator class.
	*
	* @see DecayIntegrator
	* @see PhaseSpaceChannel
	*
	* @author Peter Richardson
	*
	*/
	class PhaseSpaceMode: public Base {

	friend class PhaseSpaceChannel;

	public:

	/** @name Standard constructors and destructors. */
	//@{
	/**
	* The default constructor.
	*/
	PhaseSpaceMode() : maxWeight_(0.), partial_(-1),
	testOnShell_(false), eMax_(ZERO),
	eps_(ZERO), nRand_(0)
	{};

	/**
	* The default constructor.
	*/
	PhaseSpaceMode(tPDPtr in1, tPDVector out,
	double wMax, tPDPtr in2=tPDPtr(),
	Energy eMax=ZERO) : incoming_(make_pair(in1,in2)),
	maxWeight_(wMax),
	outgoing_(out), partial_(-1),
	testOnShell_(false), eMax_(eMax),
	eps_(ZERO), nRand_(0)
	{};
	//@}

	public:

	/**
	* Generate the decay.
	* @param intermediates Whether or not to generate the intermediate particle
	* in the decay channel.
	* @param cc Whether we are generating the mode specified or the charge
	* conjugate mode.
	* @param inpart The incoming particle.
	* @return The outgoing particles.
	*/
	ParticleVector generateDecay(const Particle & inpart,
	tcDecayIntegrator2Ptr decayer,
	bool intermediates,bool cc);

	/**
	* Add a new channel.
	* @param channel A pointer to the new PhaseChannel
	*/
	void addChannel(PhaseSpaceChannel channel) {
	channel.init(this);
	channels_.push_back(channel);
	}

	/**
	* Reset the properties of one of the intermediate particles. Only a specific channel
	* is reset.
	* @param ichan The channel to reset.
	* @param part The ParticleData object of the particle to reset
	* @param mass The mass of the intermediate.
	* @param width The width of gthe intermediate.
	*/
	void resetIntermediate(int ichan, tcPDPtr part, Energy mass, Energy width) {
	if(!part) return;
	channels_[ichan].resetIntermediate(part,mass,width);
	}

	/**
	* Reset the properties of one of the intermediate particles. All the channels
	* are reset.
	* @param part The ParticleData object of the particle to reset
	* @param mass The mass of the intermediate.
	* @param width The width of gthe intermediate.
	*/
	void resetIntermediate(tcPDPtr part, Energy mass, Energy width) {
	for(PhaseSpaceChannel & channel : channels_)
	channel.resetIntermediate(part,mass,width);
	}

	/**
	* The phase-space channels
	*/
	const vector<PhaseSpaceChannel> & channels() const {return channels_;}

	/**
	* Set the weights
	*/
	void setWeights(const vector<double> & wgts) {
	assert(wgts.size()==channels_.size());
	for(unsigned int ix=0;ix<channels_.size();++ix)
	channels_[ix].weight(wgts[ix]);
	}

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

	public :

	/**
	* Initialize the phase space
	*/
	void init();

	/**
	* Initialisation before the run stage
	*/
	void initrun();

	/**
	* Get the maximum weight for the decay.
	* @return The maximum weight.
	*/
	double maxWeight() const {return maxWeight_;}

	/**
	* Set the maximum weight for the decay.
	* @return The maximum weight.
	*/
	void maxWeight(double wgt) const {maxWeight_=wgt;}

	/**
	* Initialise the phase space.
	* @param init Perform the initialization.
	*/
	Energy initializePhaseSpace(bool init, tcDecayIntegrator2Ptr decayer,
	bool onShell=false);

	/**
	* The incoming particles
	*/
	pair<PDPtr,PDPtr> incoming() const {return incoming_;}

	/**
	* Access to the outging particles.
	* @return A pointer to the ParticleData object.
	*/
	tPDVector outgoing() const {return outgoing_;}

	/**
	* Number of outgoing particles.
	* @return The number of outgoing particles.
	*/
	unsigned int numberOfParticles() const {return outgoing_.size();}

	/**
	* Set the partial width to use for normalization. This is the partial width
	* in the WidthGenerator object.
	* @param in The partial width to use.
	*/
	void setPartialWidth(int in) {partial_=in;}

	/**
	* Access to the epsilon parameter
	*/
	Energy epsilonPS() const {return eps_;}

	public :

	/**
	* A friend operator to allow the mode to be outputted for debugging purposes.
	*/
	friend ostream & operator<<(ostream & os, const PhaseSpaceMode & mode) {
	os << "The mode has " << mode.channels_.size() << " channels\n";
	if(mode.incoming_.second==PDPtr())
	os << "This is a mode for the decay of " << mode.incoming_.first->PDGName() << " to ";
	else
	os << "This is a mode for " << mode.incoming_.first->PDGName() << ", "
	<< mode.incoming_.second->PDGName() << " to ";
	for(tPDPtr out : mode.outgoing_) os << out->PDGName() << " ";
	os << "\n";
	for(const PhaseSpaceChannel & channel : mode.channels_) os << channel;
	return os;
	}

	private:

	/**
	* Return the weight for a given phase-space point.
	- * @param ichan The channel to generate the weight for.
	- * @param in The incoming particle.
	- * @param particles The outgoing particles.
	- * @param first Whether or not this is the first call for initialisation
	+ * @param in The momentum of the incoming particle
	+ * @param momenta The momenta of the outgoing particles
	+ * @param onShell Whether or not to force the intermediates to be on-shell
	* @return The weight.
	*/
	- Energy weight(int & ichan, const Particle & in,
	+ Energy weight(int & ichan, const Lorentz5Momentum & in,
	vector<Lorentz5Momentum> & momenta,
	bool onShell=false) const {
	ichan=0;
	// flat phase-space
	if(channels_.empty())
	return flatPhaseSpace(in,momenta,onShell);
	// multi-channel
	else
	return channelPhaseSpace(ichan,in,momenta,onShell);
	}

	/**
	* Return the weight and momenta for a flat phase-space decay.
	- * @param inpart The incoming particle.
	- * @param outpart The outgoing particles.
	+ * @param in The momentum of the incoming particle
	+ * @param momenta The momenta of the outgoing particles
	+ * @param onShell Whether or not to force the intermediates to be on-shell
	* @return The weight.
	*/
	- Energy flatPhaseSpace(const Particle & inpart,
	+ Energy flatPhaseSpace(const Lorentz5Momentum & in,
	vector<Lorentz5Momentum> & momenta,
	bool onShell=false) const;

	/**
	* Generate a phase-space point using multichannel phase space.
	- * @param in The incoming particle.
	- * @param particles The outgoing particles.
	+ * @param ichan The channel to use
	+ * @param in The momentum of the incoming particle
	+ * @param momenta The momenta of the outgoing particles
	+ * @param onShell Whether or not to force the intermediates to be on-shell
	* @return The weight.
	*/
	- Energy channelPhaseSpace(int & ichan, const Particle & in,
	+ Energy channelPhaseSpace(int & ichan, const Lorentz5Momentum & in,
	vector<Lorentz5Momentum> & momenta,
	bool onShell=false) const;

	/**
	* Generate the masses of the external particles.
	* @param inmass The mass of the decaying particle.
	* @param wgt The weight for the masses.
	* @return The masses.
	*/
	vector<Energy> externalMasses(Energy inmass,double & wgt, bool onShell) const;

	/**
	* Construct the vertex for spin corrections
	* @param in The incoming particle.
	* @param out The outgoing particles.
	*/
	void constructVertex(const Particle & in, const ParticleVector & out,
	tcDecayIntegrator2Ptr decayer) const;

	/**
	* Fill the stack
	*/
	void fillStack() {
	assert(rStack_.empty());
	for(unsigned int ix=0;ix<nRand_;++ix) rStack_.push(UseRandom::rnd());
	}

	private:

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

	private:

	/**
	* The incoming particles
	*/
	pair<PDPtr,PDPtr> incoming_;

	/**
	* The phase-space channels
	*/
	vector<PhaseSpaceChannel> channels_;

	/**
	* The maximum weight
	*/
	mutable double maxWeight_;

	/**
	* The external particles
	*/
	tPDVector outgoing_;
	tPDVector outgoingCC_;

	/**
	* Which of the partial widths of the incoming particle to use
	*/
	int partial_;

	/**
	* The width generator for the incoming particle.
	*/
	cGenericWidthGeneratorPtr widthGen_;

	/**
	* The mass generators for the outgoing particles.
	*/
	vector<cGenericMassGeneratorPtr> massGen_;

	/**
	* Whether to check on-shell or off-shell kinematics
	* in doinit, if on-shell off-shell is tested in initrun
	*/
	bool testOnShell_;

	/**
	* The maximum energy for the mode
	*/
	Energy eMax_;

	/**
	* Cut-off
	*/
	Energy eps_;

	/**
	* Number of random numbers required
	*/
	unsigned int nRand_;

	/**
	* Stack for the random numbers
	*/
	mutable stack<double> rStack_;
	};

	}

	#endif /* Herwig_PhaseSpaceMode_H */

File Metadata

Mime Type: text/x-diff
Expires: Mon, Jan 20, 9:13 PM (21 h, 54 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4237765
Default Alt Text: (28 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions