External Decays

DecayHandler is a base class for the external handling of decays. It is intended for normal particle decays, primarily B mesons and tau, and cannot be used to redirect decays of heavy resonances like t or Z^0. The user-written derived class is called if a pointer to it has been given with the pythia.decayPtr() method, where it also is specified which particles it will be called for. This particle information is accessible with the doExternalDecay() method.

The main pure virtual method in DecayHandler to do the decay is:

virtual bool DecayHandler::decay(vector<int>& idProd, vector<double>& mProd, vector<Vec4>& pProd, int iDec, const Event& event)  
argument idProd : is a list of particle PDG identity codes,
argument mProd : is a list of their respective masses (in GeV), and
argument pProd : is a list of their respective four-momenta.

At input, these vectors each have size one, so that idProd[0], mProd[0] and pProd[0] contain information on the particle that is to be decayed. At output, the vectors should have increased by the addition of all the decay products. Even if initially defined in the rest frame of the mother, the products should have been boosted so that their four-momenta add up to the pProd[0] of the decaying particle.

Should it be of interest to know the prehistory of the decaying particle, e.g. to set some helicity information affecting the decay angular distribution, the full event record is available read-only, with info in which slot iDec the decaying particle is stored.

The routine should return true if it managed the decay and false otherwise, in which case Pythia will try to do the decay itself. This e.g. means you can choose to do some decay channels yourself, and leave others to Pythia. To avoid double-counting, the channels you want to handle should be switched off in the Pythia particle database. In the beginning of the external decay method you should then return false with a probability given by the sum of the branching ratios for those channels you do not want to handle yourself.

Note that the decay vertex is always set by Pythia, and that B-Bbar oscillations have already been taken into account, if they were switched on. Thus idProd[0] may be the opposite of event[iDec].id(), where the latter provides the code at production.

One limitation of the method above is that it is only intended for one decay step, not for a sequential decay chain. (At least not for displaying such intermediate steps.) That is, the control for any subsequent decays returns to PYTHIA. If you want to avoid this another method exists, with one extra argument:

virtual bool DecayHandler::chainDecay( vector<int>& idProd, vector<int>& motherProd, vector<double>& mProd, vector<Vec4>& pProd, int iDec, const Event& event)  
argument motherProd : is a list of the indices of the mother,
and the other arguments are as above.

Here the new motherProd vector also has size one at input, with motherProd[0] = 0. At output it should have increaed in size in the same way as the other arrays. Particles that come directly from the mother should have value 0, whereas secondary decay products should have the index of the mother in the arrays. To simplify parsing, particles having the same mother should be placed consecutively in the arrays, and daughters can not be put before their mothers. When the particles are transferred to the standard event record, the full mother-daughter relations will be reconstructed from the new array, and any particle with daughters will be considered to have decayed. For long-lived intermediate particles also vertex information will take this into account. User-selected secondary decay channels will be accepted as they are, however, without any knowledge whether the user has allowed for particle-antiparticle oscillations before that decay. Therefore a simple exponential decay time will be used to find secondary vertices.

While primarily intended for sequential decays, of course the chainDecay method can be used also for simple decays in one step, and is then equivalent with decay one. This is useful if a particle species has some decay channels that lead to sequential decays whereas others do not. During code execution it is first checked whether chainDecay can do the decay, and if not decay is offered to. By default chainDecay returns false, so if you only overload decay it will be called. If you want to you can choose to handle the decays of some particles in one of the methods and other particles in the other method, so long as you return false for those decays you do not handle.

The choice of which method to use can be done by the user in the method

virtual bool doChainDecay(int idMother)  
if true then chainDecay is called, else decay. This method is called for each particle that is to be decayed externally.
argument idMother : is the PDG identity code of the decaying mother.

A sample test program is available in main17.cc, providing a simple example of how to use this facility.


The external B and C-hadron decay program EvtGen performs a chain of decays, rather than single particle decays, to propagate helicity information throughout the chain. Consequently, EvtGen cannot be simply interfaced via the DecayHandler class. A special class, EvtGenDecays, is provided in Pythia8Plugins which can be called after an event has been generated, to perform all remaining decays via EvtGen. An example of how to use this class is provided in main48.cc. A more detailed discussion of some physics considerations, notably event weights for forced decays, can be found in this note.