Implement New Showers

In case you want to replace the PYTHIA initial- and final-state showers by your own, it is possible but not trivial. The point is that multiple interactions (MI), initial-state radiation (ISR) and final-state radiation (FSR) in general appear in one single interleaved sequence of decreasing pT values. Therefore shower replacements would have to be able to play the game by such rules, as we will outline further below. Of course, this still leaves the field open exactly how to define what to mean by pT, how to handle recoil effects, how the colour flow is affected, and so on, so there is certainly room for alternative showers.

For the moment we assume you want to keep the MI part of the story unchanged, and make use of the existing beam-remnants (BR) machinery. If you want to replace both MI, ISR, FSR and BR then you had better replace the whole PartonLevel module of the code. If, in addition, you want to produce your own hard processes, then you only need the hadron-level standalone part of the machinery.

In order to write replacement codes for ISR and/or FSR it is useful to be aware of which information has to be shared between the different components, and which input/output structure is required of the relevant methods. For details, nothing beats studying the existing code. However, here we provide an overview, that should serve as a useful introduction.

It should be noted that we here primarily address the problem in its full generality, with interleaved MI, ISR and FSR. There exists an option TimeShower:interleave = off where only MI and ISR would be interleaved and FSR be considered after these two, but still before BR. Most of the aspects described here would apply also for that case. By contrast, resonance decays are only considered after all the four above components, and timelike showers in those decays would never be interleaved with anything else, so are much simpler to administrate.

Therefore the pythia.setShowerPtr( timesDecPtr, timesPtr, spacePtr) method allows two separate pointers to be set to instances of derived TimeShower classes. The first is only required to handle decays, say of Z^0 or Upsilon, with no dependence on beam remnants or ISR. The second, as well as spacePtr, has to handle the interleaved evolution of MI, ISR and FSR. Therefore you are free to implement only the first, and let the PYTHIA default showers take care of the latter two. But, if you wanted to, you could also set timesDecPtr = 0 and only provide a timesPtr, or only a spacePtr. If your timelike shower does both cases, the first two pointers can agree. The only tiny point to take into account then is that init( beamAPtr, beamBPtr) is called twice, a first time to timesDecPtr with beam pointers 0, and a second time to timesPtr with nonvanishing beam pointers.

The event record

Obviously the main place for sharing information is the event record, specifically the Event event member of Pythia, passed around as a reference. It is assumed you already studied how it works, so here we only draw attention to a few aspects of special relevance.

One basic principle is that existing partons should not be overwritten. Instead new partons should be created, even when a parton only receives a slightly shifted momentum and for the rest stays the same. Such "carbon copies" by final-state branchings should be denoted by both daughter indices of the original parton pointing to the copy, and both mother indices of the copy to the original. If the copy instead is intended to represent an earlier step, e.g. in ISR backwards evolution, the role of mothers and daughters is interchanged. The event.copy( iCopy, newStatus) routine can take care of this tedious task; the sign of newStatus tells the program which case to assume.

To make the event record legible it is essential that the status codes are selected appropriately to represent the reason why each new parton is added to the record. Also remember to change the status of a parton to be negative whenever an existing parton is replaced by a set of new daughter partons.

Another important parton property is scale(), which does not appear in the normal event listing, but only if you use the extended Event:listScaleAndVertex = on option. This property is supposed to represent the production scale (in GeV) of a parton. In the current FSR and ISR algorithms it is used to restrict from above the allowed pT values for branchings of this particular parton. Beam remnants and other partons that should not radiate are assigned scale 0.

The subsystems

One aspect that complicates administration is that an event can contain several subsystems, consisting of one MI and its associated ISR and FSR, which to first approximation are assumed to evolve independently, but to second are connected by the interleaved evolution. The partons of a given subsystem therefore do not have to be stored consecutively. Some simple vectors put inside the event record can be used to keep track of the current position, i.e. index iPos for parton event[iPos], of all partons of all systems.

The number of systems is given by sizeSystems(), with systems numbered 0 through sizeSystems() - 1. The size of a given subsystem iSys is given by sizeSystem(iSys), with analogous numbering. The method getInSystem( iSys, iMem) returns the position iPos of the iMem'th member of the iSys'th system. For each system, the slots iMem = 0 and 1 are intended for the incoming partons from beam A and B, respectively. If there are no beams, such as in resonance decays, these should be assigned value 0. Slots 2 onwards are intended for all the outgoing partons. These latter partons are not guaranteed to appear in any particular order.

New systems are created from the hard process and by the MI, not from any of the other components, by the newSystem() method, which returns the iSys code of the new system. Both FSR and ISR modify the position of partons, however. Since an FSR or ISR branching typically implies a new state with one more parton than before, the method addToSystem( iSys, iPos) appends a new slot to the given system and stores the iPos value there. Furthermore, in a branching, several existing partons may also be moved to new slots. The method setInSystem( iSys, iMem, iPos) replaces the current iPos value in the given slot by the provided one. If you do not know the iMem value, replaceInSystem( iSys, iPosOld, iPosNew) will replace any ocurence of iPosOld by iPosNew wherever it is found in the iSys'th system. In a FSR 1 -> 2 branching it is irrelevant which parton position you let overwrite the existing slot and which is added to the end of the system.

Finally, clearSystems() empties the information, and listSystems() provides a simple listing of all the iPos values by system, intended for debugging purposes only.

The system information must be kept up-to-date. The MI component only writes, but both ISR, FSR and BR make extensive use of the existing information. As an example, the introduction of primordial kT in the beam remnants will fail if the information on which final-state partons belong to which system is out-of-date.

Currently the system information is kept throughout the continued history of the event. Resonance decays create new systems, appended to the existing ones. This could be useful during the hadronization stage, to collect the partons that belong to a resonace with preserved mass when a small string collapses to one particle, but is not yet used for that.

The beams

The different subsystems are tied together by them sharing the same initial beam particles, and thereby being restricted by energy-momentum and flavour conservation issues. The information stored in the two beam particles, here called beamA and beamB, is therefore as crucial to keep correct as the above subsystem list.

Both beam objects are of the BeamParticle class. Each such object contains a vector with the partons extracted from it. The number of such partons, beamX.size() (X = A or B), of course is the same as the above number of subsystems in the event record. (The two diverge at the BR step, where further beam remnants are added to the beams without corresponding to new subsystems.) The individual partons are accessed by an overloaded indexing operator to a vector of ResolvedParton objects. The iPos() property corresponds to the iPos one above, i.e. providing the position in the main event record of a parton. In particular, beamA[iSys].iPos() = event.getInSystem( iSys, 0) and beamB[iSys].iPos() = event.getInSystem( iSys, 1). Whereas thus the indices of the two incoming partons to a subsystem are stored in two places, the ones of the outgoing partons only appear in the system part of the Event class.

Just as the subsystems in Event must be updated, so must the information in the two BeamParticle's, e.g. with methods beamX[iSys].iPos( iPosIn) when an incoming parton is replaced by a new one in line iPosIn. Furthermore the new parton identity should be set by beamX[iSys].id( idIn) and the new x energy-momentum fraction by beamX[iSys].x( xIn). The three can be combined in one go by beamX[iSys].update( iPosIn, idIn, xIn).

To be specific, it is assumed that, at each step, the two incoming partons are moving along the +-z axis and are massless. Since the event is constructed in the c.m. frame of the incoming beams this implies that x = 2 E / E_cm. If the x values are not defined accordingly or not kept up-to-date the BR treatment will not conserve energy-momentum.

In return, the BeamParticle objects give access to some useful methods. The beamX.xf( id, x, Q2) returns the standard PDF weight x f_id(x, Q^2). More intererstingly, beamX.xfISR( iSys, id, x, Q2) returns the modified weight required when several subsystems have to share the energy and flavours. Thus iSys is added as an extra argument, and the momentum already assigned to the other subsystems is not available for evolution, i.e. the maximal x is correspondingly smaller than unity. Also flavour issues are handled in a similar spirit.

An additional complication is that a parton can be either valence or sea, and in the latter case the BR treatment also distinguishes companion quarks, i.e. quark-antiquark pairs that are assumed to come from the same original g -> q qbar branching, whether perturbative or not. This can be queried either with the beamX[iSys].companion() method, for detailed information, or with the beamX[iSys].isValence(), beamX[iSys].isUnmatched() and beamX[iSys].isCompanion() metods for yes/no answers whether a parton is valence, unmatched sea or matched sea. This choice should affect the ISR evolution; e.g., a valence quark cannot be constructed back to come from a gluon.

To keep this info up-to-date, the beamX.pickValSeaComp() method should be called whenever a parton of a new flavour has been picked in the ISR backwards evolution, but not if the flavour has not been changed, since then one should not be allowed to switch back and forth between the same quark being considered as valence or as sea. Since the pickValSeaComp() method makes use of the current parton-density values, it should be preceded by a call to beamX.xfISR( iSys, id, x, Q2), where the values in the call are the now finally accepted ones for the newly-found mother. (Such a call is likely to have been made before, during the evolution, but is not likely to be the most recent one, i.e. still in memory, and therefore had better be redone.)

Most of the issues in this section are related to the ISR algorithm, i.e. it is possible to write an FSR algorithm that is completely decoupled. Alternatively the FSR may change the position where an incoming parton is stored, or its assumed momentum, e.g. by recoil effects inside dipoles stretched from the scattered back to the incoming partons. In that case some of the methods above would have to be used also inside the FSR algorithm.

The TimeShower interface

If you want to replace the TimeShower class this would involve replacing the following methods.

method TimeShower()
The constructor does not need to do anything.

method void init( BeamParticle* beamAPtrIn = 0, BeamParticle* beamBPtrIn = 0)
You have to store your local copy of the pointers to these objects, since they have to be used during the generation, as explained above. The pointers could be zero; e.g. a local copy of TimeShower is created to handle showers in decays such as Upsilon -> q qbar from inside the ParticleDecays class. This is also the place to do initialization of whatever parameters you plan to use, e.g. by reading in them from a user-accessible database like the Settings one.

method double enhancePTmax()
Relative to the default pT_max evolution scale of the process, it may still be convenient to vary the matching slightly for the hardest interaction in an event, to probe the sensitivity to such details. The base-class implementation returns the value of the TimeShower:pTmaxFudge parameter.

method int shower( int iBeg, int iEnd, Event& event, double pTmax)
This is an all-in-one call for shower evolution, and as such cannot be used for the normal interleaved evolution, where only the routines below are used. It also cannot be used in resonance decays that form part of the hard process, since there the user hooks insert a potential veto step. Currently this routine is therefore only used in the hadron-level decays, e.g. Upsilon -> g g g.
iBeg and iEnd is the position of the first and last parton of a separate system, typically produced by a resonance decay. Such a system only evolves in isolation, and in particular does not relate to the beams.
The pTmax value sets the maximum scale for evolution, but normally you would restrict that further for each individual parton based on its respective scale value.
The routine is expected to return the number of FSR branchings that were generated, but only for non-critical statistics purposes.
Since the real action typically is delegated to the routines below, it may well be that the existing code need not be replaced.

method void prepare( int iSys, Event& event)
This method is called immediately after a new interaction (or the products of a resonance decay) has been added, and should then be used to prepare the subsystem of partons for subsequent evolution. In the current code this involves identifying all colour and charge dipole ends: the position of radiating and recoiling partons, maximum pT scales, possible higher-order matrix elements matchings to apply, and so on.
The iSys parameter specifies which parton system is to be prepared. It is used to extract the set of partons to be treated, with rules as described in the above section on subsystems. Specifically, the first two partons represent the incoming state, or are 0 for resonance decays unrelated to the beams, while the rest are not required to be in any particular order.

method void update( int iSys, Event& event)
This method is called immediately after a spacelike branching in the iSys'th subsystem. Thus the information for that system is out-of-date, while that of the others is unchanged. If you want, you are free to throw away all information for the affected subsystem and call prepare( iSys, event) to create new one. Alternatively you may choose only to update the information that has changed.

method double pTnext( Event& event, double pTbegAll, double pTendAll)
This is the main driver routine for the downwards evolution. A new pT is to be selected based on the current information set up by the routines above, and along with that a branching parton or dipole. The pTbegAll scale is the maximum scale allowed, from which the downwards evolution should be begun (usually respecting the maximum scale of each individual parton). If no emission is found above pTendAll (and above the respective shower cutoff scales) then 0. should be returned and no emissions will be allowed. Both scales can vary from one event to the next: if a scale has already been selected for MI or ISR it makes no sense to look for a scale smaller than that from FSR, since it would not be able to compete, so pTendAll is set correspondingly. As it happens, FSR is tried before ISR and MI in the interleaved evolution, but this is an implementational detail that could well change.
Typically the implementation of this routine would be to set up a loop over all possible radiating objects (dipoles, dipole ends, ...), for each pick its possible branching scale and then pick the one with largest scale as possible winner. At this stage no branching should actually be carried out, since MI, ISR and FSR still have to be compared to assign the winner.

method bool branch( Event& event)
This method will be called once FSR has won the competition with MI and ISR to do the next branching. The candidate branching found in the previous step should here be carried out in full. The pre-branching partons should get a negative status code and new replacement ones added to the end of the event record. Also the subsystem information should be updated, and possibly also the beams.
Should some problem be encountered in this procedure, e.g. if some not-previously-considered kinematics requirement fails, it is allowed to return false to indicate that no branching could be carried out.

method int system()
This method is not virtual. If a branching is constructed by the previous routine this tiny method should be able to return the number of the selected subsystem iSysSel where it occured, so that the spacelike shower can be told which system to update, if necessary. Therefore iSysSel must be set in branch (or already in pTnext).

method void list( ostream& os = cout)
This method is not at all required. In the current implementation it outputs a list of all the dipole ends, with information on the respective dipole. The routine is not called anywhere in the public code, but has been inserted at various places during the development/debug phase.

The SpaceShower interface

If you want to replace the SpaceShower class this would involve replacing the following methods. You will find that much of the story reminds of TimeShower above, and actually some cut-and-paste of text is involved. In some respects the description is simpler, since there are no special cases for resonance decays and non-interleaved evolution. Thus there is no correspondence to the TimeShower::shower(...) routine.

method SpaceShower()
The constructor does not need to do anything.

method void init( BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn)
You have to store your local copy of the pointers to these objects, since they have to be used during the generation, as explained above. This is also the place to do initialization of whatever parameters you plan to use, e.g. by reading in them from a user-accessible database like the Settings one.

method bool limitPTmax( Event& event, double Q2Fac = 0., double Q2Ren = 0.)
The question is whether the ISR should be allowed to occur at larger scales than the hard process it surrounds. This is process-dependent. For instance, if the hard process is Z^0 production we know that ISR should be allowed to go right up to the kinematical limit. If it is a 2 -> 2 QCD process the ISR should not exceed the scale of the hard process, since if so one would doublecount. The SpaceShower:pTmaxMatch switch allows you to force the behaviour, or else to program your own logic. The current implementation limits pT whenever the final state contains a quark (except top), gluon or photon, since then the danger of doublecounting is there. You may replace by your own logic, or leave as is.
The internal PYTHIA implementation also allows intermediate options, where emissions can go up to the kinematical limit but be dampened above the factorization or renormalization scale. Therefore the (square of the) latter two are provided as optional input parameters.

method double enhancePTmax()
When the above method limits pT_max to the scale of the process, it may still be convenient to vary the matching slightly for the hardest interaction in an event, to probe the sensitivity to such details. The base-class implementation returns the value of the SpaceShower:pTmaxFudge parameter.

method void prepare( int iSys, Event& event, bool limitPTmax = true)
This method is called immediately after a new interaction has been added, and should then be used to prepare the subsystem of partons for subsequent evolution. In the current code this involves identifying the colour and charge dipole ends: the position of radiating and recoiling partons, maximum pT scales, and possible higher-order matrix elements matchings to apply. Depending on what you have in mind you may choose to store slightly different quantities. You have to use the subsystem information described above to find the positions of the two incoming partons (and the outgoing ones) of the system, and from there the scales at which they were produced.
The limitPTmax input agrees with the output of the previous method for the hardest process, and is always true for subsequent MI, since there an unlimited pT for sure would lead to doublecounting.

method void update( int iSys, Event& event)
This method is called immediately after a timelike branching in the iSys'th subsystem. Thus the information for that system may be out-of-date, and to be updated. For the standard PYTHIA showers this routine does not need to do anything, but that may be different in another implementation.

method double pTnext(( Event& event, double pTbegAll, double pTendAll, int nRadIn = -1)
This is the main driver routine for the downwards evolution. A new pT is to be selected based on the current information set up by the routines above, and along with that a branching parton or dipole. The pTbegAll scale is the maximum scale allowed, from which the downwards evolution should be begun (usually respecting the maximum scale of each individual parton). If no emission is found above pTendAll (and above the respective shower cutoff scales) then 0. should be returned and no emissions will be allowed. Both scales can vary from one event to the next: if a scale has already been selected for MI or ISR it makes no sense to look for a scale smaller than that from FSR, since it would not be able to compete, so pTendAll is set correspondingly. As it happens, FSR is tried before ISR and MI in the interleaved evolution, but this is an implementational detail that could well change.
Typically the implementation of this routine would be to set up a loop over all possible radiating objects (dipoles, dipole ends, ...), for each pick its possible branching scale and then pick the one with largest scale as possible winner. At this stage no branching should actually be carried out, since MI, ISR and FSR still have to be compared to assign the winner.
The final input nRadIn provides the total number of ISR and FSR emissions already generated in the event, and so allows a special treatment for the very first emission, if desired.