Page MenuHomeHEPForge
Files F8724006

No OneTemporary
Actions

View Options

diff --git a/Decay/PhaseSpaceMode.cc b/Decay/PhaseSpaceMode.cc
--- a/Decay/PhaseSpaceMode.cc
+++ b/Decay/PhaseSpaceMode.cc
@@ -1,427 +1,427 @@
// -*- 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);
}
void PhaseSpaceMode::persistentInput(PersistentIStream & is, int) {
is >> incoming_ >> channels_ >> maxWeight_ >> outgoing_ >> outgoingCC_
>> partial_ >> widthGen_ >> massGen_ >> testOnShell_
>> iunit(eMax_,GeV);
}
// 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());
// 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
wgt = pre*weight(ichan,inrest,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) {
- decayer->me2(ix,inrest,!cc ? outgoing_ : outgoingCC_,
- momenta,DecayIntegrator2::Calculate);
+ mewgts[ix]=decayer->me2(ix,inrest,!cc ? outgoing_ : outgoingCC_,
+ momenta,DecayIntegrator2::Calculate);
total+=mewgts[ix];
}
// randomly pick a channel
total *= UseRandom::rnd();
- iChannel_ = 0;
+ int iChannel = -1;
do {
- total-=mewgts[iChannel_];
- ++iChannel_;
+ ++iChannel;
+ total-=mewgts[iChannel];
}
- while(iChannel_<channels_.size() && total>0.);
+ while(iChannel<int(channels_.size()) && total>0.);
// apply boost
for(tPPtr part : output) part->transform(boostFromRest);
- channels_[iChannel_].generateIntermediates(cc,inpart,output);
+ 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,
vector<Lorentz5Momentum> & momenta,
bool onShell) const {
double wgt(1.);
// masses of the particles
vector<Energy> mass = externalMasses(inpart.mass(),wgt,onShell);
// two body decay
double ctheta,phi;
Kinematics::generateAngles(ctheta,phi);
if(! Kinematics::twoBodyDecay(inpart.momentum(), 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();
}
// 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);
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());
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
wgt = pre*weight(dummy,*inpart,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 {
wgt = pre*weight(ichan,*inpart,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,
vector<Lorentz5Momentum> & momenta,
bool onShell) const {
double wgt(UseRandom::rnd());
// 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);
// 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;
}
diff --git a/Decay/PhaseSpaceMode.h b/Decay/PhaseSpaceMode.h
--- a/Decay/PhaseSpaceMode.h
+++ b/Decay/PhaseSpaceMode.h
@@ -1,369 +1,370 @@
// -*- 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), iChannel_(999), eMax_(ZERO)
+ testOnShell_(false), eMax_(ZERO)
{};
/**
* 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), iChannel_(999), eMax_(eMax)
+ testOnShell_(false), eMax_(eMax)
{};
//@}
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() {
// 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();
}
void 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();
}
}
/**
* 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> const incoming() {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;}
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
* @return The weight.
*/
Energy weight(int & ichan, const Particle & 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.
* @return The weight.
*/
Energy flatPhaseSpace(const Particle & inpart,
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.
* @return The weight.
*/
Energy channelPhaseSpace(int & ichan, const Particle & 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;
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 selected channel
- */
- mutable unsigned int iChannel_;
-
- /**
* The maximum energy for the mode
*/
Energy eMax_;
};
}
#endif /* Herwig_PhaseSpaceMode_H */

File Metadata

Mime Type
text/x-diff
Expires
Mon, Jan 20, 10:21 PM (1 d, 7 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4242604
Default Alt Text
(26 KB)

Event Timeline