Page MenuHomeHEPForge
Files F7878610

No OneTemporary
Actions

View Options

diff --git a/Decay/DecayIntegrator2.h b/Decay/DecayIntegrator2.h
--- a/Decay/DecayIntegrator2.h
+++ b/Decay/DecayIntegrator2.h
@@ -1,477 +1,477 @@
// -*- C++ -*-
#ifndef Herwig_DecayIntegrator2_H
#define Herwig_DecayIntegrator2_H
//
// This is the declaration of the DecayIntegrator2 class.
//
#include "DecayIntegrator2.fh"
#include "HwDecayerBase.h"
#include "PhaseSpaceMode.fh"
#include "Herwig/PDT/WidthCalculatorBase.fh"
#include "Radiation/DecayRadiationGenerator.h"
#include <Herwig/Decay/DecayVertex.h>
namespace Herwig {
using namespace ThePEG;
/** \ingroup Decay
* \class DecayIntegrator2
* \brief Main class for Decayers implementing multi-channel phase space integration.
* \author Peter Richardson
*
* This class is designed to be the base class for Herwig decays including
* the implementation of a multichannel decayer or n-body phase space decays.
*
* The <code>DecayIntegrator2</code> class inherits from ThePEG's Decayer class
* and makes use of the <code>PhaseSpaceMode</code> class to specify a number
* of decay modes.
*
* Additional modes can be added using the addMode method. In practice the
* phase space channels for a particular mode are usually constructed in the
* doinit member of a Decayer and then the modes added to the Decayer.
*
* For the majority of the decays currently implemented the
* phase-space integration has been optimised and the maximum weight set.
* If the parameters of the decay model are changed the Initialize interface
* can be used to optimise the integration and calculate the maximum weight.
*
* In classes inheriting from this the me2() member which gives the matrix element
* squared must be implemented. This should be combined with the setting of the
* phase space channels, and the setting of which channels to use and their
* initial weights in the doinit() member. The different decay modes should then
* be initialized in the initrun() member if needed. The generate member can then
* be called from the decay() member to generate a phase-space configuration for a
* decay.
*
* @see DecayPhaseSpaceMode
* @see DecayPhaseSpaceChannel
* @see \ref DecayIntegrator2Interfaces "The interfaces"
* defined for DecayIntegrator2.
*/
class DecayIntegrator2: public HwDecayerBase {
public:
/**
* and DecayPhaseMode
*/
friend class PhaseSpaceMode;
/**
* Enum for the matrix element option
*/
enum MEOption {Initialize,Calculate,Terminate};
public:
/**
* The default constructor.
*/
DecayIntegrator2() : nIter_(10), nPoint_(10000), nTry_(500),
generateInter_(false), iMode_(-1),
realME_(false), virtualME_(false), eps_(ZERO)
{}
public:
/**
* Check if this decayer can perfom the decay for a particular mode.
* Uses the modeNumber member but can be overridden
* @param parent The decaying particle
* @param children The decay products
*/
virtual bool accept(tcPDPtr parent, const tPDVector & children) const {
bool cc;
return modeNumber(cc,parent,children)>=0;
}
/**
* For a given decay mode and a given particle instance, perform the
* decay and return the decay products. As this is the base class this
* is not implemented.
* @return The vector of particles produced in the decay.
*/
virtual ParticleVector decay(const Particle & parent,
const tPDVector & children) const;
/**
* Which of the possible decays is required
* @param cc Is this mode the charge conjugate
* @param parent The decaying particle
* @param children The decay products
*/
virtual int modeNumber(bool & cc, tcPDPtr parent,
const tPDVector & children) const = 0;
/**
* The mode being used for this decay
*/
int imode() const {return iMode_;}
/**
* Add a phase-space mode to the list
* @param mode The mode being added.
*/
void addMode(PhaseSpaceModePtr mode) const;
/**
* Return the matrix element squared for a given mode and phase-space channel.
* This function is purely virtual and must be implemented in classes inheriting
* from DecayIntegrator.
* @param ichan The channel we are calculating the matrix element for.
* @param part The decaying Particle.
* @param decay The particles produced in the decay.
* @param opt Option for the calculation of the matrix element
* @return The matrix element squared for the phase-space configuration.
*/
virtual double me2(const int ichan, const Particle & part,
const tPDVector & outgoing,
const vector<Lorentz5Momentum> & momenta,
MEOption opt) const=0;
virtual void constructSpinInfo(const Particle & part, ParticleVector outgoing) const = 0;
/**
* Output the setup information for the particle database
* @param os The stream to output the information to
* @param header Whether or not to output the information for MySQL
*/
virtual void dataBaseOutput(ofstream & os,bool header) const;
/**
* Set the code for the partial width. Finds the partial width in the
* GenericWidthGenerator class which corresponds to the decay mode.
* @param dm The DecayMode
* @param imode The mode.
*/
void setPartialWidth(const DecayMode & dm, int imode);
/**
* Specify the \f$1\to2\f$ matrix element to be used in the running width calculation.
* @param mecode The code for the matrix element as described
* in the GenericWidthGenerator class.
* @param coupling The coupling for the matrix element.
* @return True or False if this mode can be handled.
*/
virtual bool twoBodyMEcode(const DecayMode &, int & mecode,
double & coupling) const {
coupling = 1.;
mecode = -1;
return false;
}
/**
* Method to return an object to calculate the 3 (or higher body) partial width
* @param dm The DecayMode
* @return A pointer to a WidthCalculatorBase object capable of calculating the width
*/
virtual WidthCalculatorBasePtr threeBodyMEIntegrator(const DecayMode & dm) const;
/**
* The matrix element to be integrated for the three-body decays as a function
* of the invariant masses of pairs of the outgoing particles.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s3 The invariant mass squared of particles 1 and 2, \f$s_3=m^2_{12}\f$.
* @param s2 The invariant mass squared of particles 1 and 3, \f$s_2=m^2_{13}\f$.
* @param s1 The invariant mass squared of particles 2 and 3, \f$s_1=m^2_{23}\f$.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The matrix element
*/
virtual double threeBodyMatrixElement(const int imode, const Energy2 q2,
const Energy2 s3, const Energy2 s2,
const Energy2 s1, const Energy m1,
const Energy m2, const Energy m3) const;
/**
* The differential three body decay rate with one integral performed.
* @param imode The mode for which the matrix element is needed.
* @param q2 The scale, \e i.e. the mass squared of the decaying particle.
* @param s The invariant mass which still needs to be integrate over.
* @param m1 The mass of the first outgoing particle.
* @param m2 The mass of the second outgoing particle.
* @param m3 The mass of the third outgoing particle.
* @return The differential rate \f$\frac{d\Gamma}{ds}\f$
*/
virtual InvEnergy threeBodydGammads(const int imode, const Energy2 q2,
const Energy2 s,
const Energy m1, const Energy m2,
const Energy m3) const;
/**
* Finds the phase-space mode corresponding to a given decay mode
* @param dm The DecayMode
*/
int findMode(const DecayMode & dm);
public:
/**
* Members for the generation of QED radiation in the decays
*/
//@{
/**
* Use the DecayRadiationGenerator to generate photons in the decay.
* @param p The Particle instance being decayed
* @param children The decay products
* @return A particle vector containing the decay products after the generation
* of photons.
*/
ParticleVector generatePhotons(const Particle & p,ParticleVector children) {
assert(false);
//return photonGen_->generatePhotons(p,children,this);
}
/**
* check if photons can be generated in the decay
*/
bool canGeneratePhotons() {return photonGen_;}
/**
* The one-loop virtual correction.
* @param imode The mode required.
* @param part The decaying particle.
* @param products The decay products including the radiated photon.
* @return Whether the correction is implemented
*/
virtual double oneLoopVirtualME(unsigned int imode,
const Particle & part,
const ParticleVector & products);
/**
* Whether or not the one loop matrix element is implemented
*/
bool hasOneLoopME() {return virtualME_;}
/**
* The real emission matrix element
* @param imode The mode required
* @param part The decaying particle
* @param products The decay products including the radiated photon
* @param iemitter The particle which emitted the photon
* @param ctheta The cosine of the polar angle between the photon and the
* emitter
* @param stheta The sine of the polar angle between the photon and the
* emitter
* @param rot1 Rotation from rest frame to frame for real emission
* @param rot2 Rotation to place emitting particle along z
*/
virtual InvEnergy2 realEmissionME(unsigned int imode,
const Particle & part,
ParticleVector & products,
unsigned int iemitter,
double ctheta, double stheta,
const LorentzRotation & rot1,
const LorentzRotation & rot2);
/**
* Whether or not the real emission matrix element is implemented
*/
bool hasRealEmissionME() {return realME_;}
//@}
public:
/**
* The output operator is a friend, this is mainly for debugging
*/
friend ostream & operator<<(ostream & os, const DecayIntegrator2 & decay);
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 Standard Interfaced functions. */
//@{
/**
* Initialize this object. Called in the run phase just before
* a run begins.
*/
virtual void doinitrun();
//@}
protected:
/**
* Set the mode being use for this decay.
*/
void imode(int in) { iMode_ = in;}
/**
* Set the helicity matrix element for the decay.
*/
void ME(DecayMEPtr in) const { matrixElement_ = in;}
/**
* The helicity amplitude matrix element for spin correlations.
*/
DecayMEPtr ME() const {return matrixElement_;}
/**
* Reset the properities of all intermediates.
* @param part The intermediate particle being reset.
* @param mass The mass of the particle.
* @param width The width of the particle.
*/
void resetIntermediate(tcPDPtr part, Energy mass, Energy width);
protected:
/**
* Methods to set variables in inheriting classes
*/
//@{
/**
* Set whether or not the intermediates are included
*/
void generateIntermediates(bool in) {generateInter_=in;}
/**
* Whether or not the one loop matrix element is implemented
*/
void hasOneLoopME(bool in) {virtualME_=in;}
/**
* Whether or not the real emission matrix element is implemented
*/
void hasRealEmissionME(bool in) {realME_=in;}
/**
* Set the epsilon parameter
*/
void epsilonPS(Energy in) {eps_=in;}
//@}
/**
* Number of decay modes
*/
unsigned int numberModes() const {return modes_.size();}
/**
* Pointer to a mode
*/
tPhaseSpaceModePtr mode(unsigned int ix) {
return modes_[ix];
}
/**
* Pointer to a mode
*/
tcPhaseSpaceModePtr mode(unsigned int ix) const {
return modes_[ix];
}
private:
/**
* The assignment operator is private and must never be called.
* In fact, it should not even be implemented.
*/
DecayIntegrator2 & operator=(const DecayIntegrator2 &);
/**
* Parameters for the integration
*/
//@{
/**
* Number of iterations for th initialization.
*/
unsigned int nIter_;
/**
* Number of points for initialisation
*/
unsigned int nPoint_;
/**
* number of attempts to generate the decay
*/
unsigned int nTry_;
/**
* List of the decay modes
*/
mutable vector<PhaseSpaceModePtr> modes_;
//@}
/**
* Whether to include the intermediates whne outputing the results.
*/
bool generateInter_;
/**
* Pointer to the object generating the QED radiation in the decay
*/
DecayRadiationGeneratorPtr photonGen_;
/**
* mode currently being generated
*/
mutable int iMode_;
/**
* The helicity matrix element for the current decay
*/
mutable DecayMEPtr matrixElement_;
/**
* Whether or not the real photon emission matrix element exists
*/
bool realME_;
/**
* Whether or not the one-loop matrix element exists
*/
bool virtualME_;
/**
* Epsilon parameter for phase-space integration
*/
Energy eps_;
-private:
+protected:
/**
* Exception for this class and those inheriting from it
*/
class DecayIntegrator2Error: public Exception {};
};
/**
* Output information on the DecayIntegrator for debugging purposes
*/
ostream & operator<<(ostream &, const DecayIntegrator2 &);
}
#endif /* Herwig_DecayIntegrator2_H */
diff --git a/Decay/PhaseSpaceMode.cc b/Decay/PhaseSpaceMode.cc
--- a/Decay/PhaseSpaceMode.cc
+++ b/Decay/PhaseSpaceMode.cc
@@ -1,424 +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);
total+=mewgts[ix];
}
// randomly pick a channel
total *= UseRandom::rnd();
iChannel_ = 0;
do {
total-=mewgts[iChannel_];
++iChannel_;
}
while(iChannel_<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,
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;
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 6:26 PM (1 d, 21 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805588
Default Alt Text
(30 KB)

Event Timeline