## 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 multiparton interactions (MPI), 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. A first example of a shower implemented within the PYTHIA context is VINCIA, which however so far only handles FSR.

For the moment we assume you want to keep the MPI part of the story unchanged, and make use of the existing beam-remnants (BR) machinery. If you want to replace both MPI, 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 MPI, ISR and FSR. There exists an option `TimeShower:interleave = off` where only MPI 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 MPI, 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 and associated information

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.

Auxiliary to the event record proper is the `PartonSystems` class, that keep track of which partons belong together in the same scattering subsystem. This information must be kept up-to-date during the shower evolution.

For initial-state showers it is also necessary to keep track of the partonic content extracted from the beams. This information is stored in the `BeamParticle` class.

### The TimeShower interface

If you want to replace the `TimeShower` class this would involve replacing the virtual methods among the following ones.

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

virtual TimeShower::~TimeShower()
The destructor does not need to do anything.

void TimeShower::initPtr(Info* infoPtr, Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtr, CoupSM* coupSMPtr, PartonSystems* partonSystemsPtr, UserHooks* userHooksPtr)
This method only imports pointers to standard facilities, and is not virtual.

virtual void TimeShower::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.

virtual bool TimeShower::limitPTmax( Event& event, double Q2Fac = 0., double Q2Ren = 0.)
The question is whether the FSR should be allowed to occur at larger scales than the hard process it surrounds. This is process-dependent, as illustrated below for the the analogous `SpaeShower::limitPTmax(...)` method, although the two kinds of radiation need not have to be modeled identically. The `TimeShower:pTmaxMatch` switch allows you to force the behaviour among three options, but you may replace by your own logic.
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.

double TimeShower::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.

virtual int TimeShower::shower( int iBeg, int iEnd, Event& event, double pTmax, int nBranchMax = 0)
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 `nBranchMax` value, if positive, gives the maximum number of allowed branchings in the call, as useful for matching studies.
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.

virtual int TimeShower::showerQED( int iBeg, int iEnd, Event& event, double pTmax)
This is a further simplified version of the `shower` method above. Currently it only handles the emission of photons in the decay of a hadron into a pair of leptons, either a charged lepton-antilepton or a lepton-neutrino pair. It is properly matched to the matrix element in the decay via a virtual photon or W^+-, respectively. It is called as part of such decays if `ParticleDecays:allowPhotonRadiation = on`, which is not the default value.

double TimeShower::pTLastInShower()
Can be used to return the pT evolution scale of the last branching in the cascade generated with the above `shower(...)` method. Is to be set in the internal `pTLastInShower` variable, and should be 0 if there were no branchings. Can be useful for matching studies.

virtual void TimeShower::prepareGlobal( Event& event)
This method resets some counters and extracts the locations of outgoing partons, in preparation of using the optional global recoil scheme. Unlike `prepare(...)` below it is only called once during the parton-level evolution, since it only relates to the hardest interaction. Is probably of no use to most people.

virtual void TimeShower::prepare( int iSys, Event& event, bool limitPTmaxIn = true, double pTfirstTrialIn = 1e9)
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.
The `limitPTmaxIn` switch conveys the choice made on maximum scale for the dipole-ends evolution associated with the hard interaction of the event (while scales for subsequent MPIs by default are set to respect pT ordering). If `true` then this scale is set by the user choice options, see the `TimeShower::limitPTmax` method above, while if `false` emissions are allowed to go up to the kinematical limit. For the two-hard-interactions scenario, the two class variables `dopTlimit1` and `dopTlimit2` instead convey the choice made.
The `pTfirstTrialIn` parameter is tentatively used in the new weak-shower part of the code. It is related to the option to shift the weak-shower evolution scale up by (some factor times) the gauge boson mass, without this backfiring on the QCD/QED evolution. Still under evolution, and normally not used.

virtual void TimeShower::rescatterUpdate( int iSys, Event& event)
This method is called immediately after rescattering in the description of multiparton interactions. Thus the information on one or several systems is out-of-date, while that of the others is unchanged. We do not provide the details here, since we presume few implementors of new showers will want to touch the technicalities involved in obtaining a description of rescattering.

virtual void TimeShower::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.

virtual double TimeShower::pTnext( Event& event, double pTbegAll, double pTendAll, bool isFirstTrialIn = false)
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 MPI 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 MPI in the interleaved evolution, but this is an implementation 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 MPI, ISR and FSR still have to be compared to assign the winner.
The `isFirstTrialIn` is `true` only for the very first emission, and can then optionally be used to switch the weak-shower evolution scale up by (some factor times) the gauge boson mass.

virtual bool TimeShower::branch( Event& event, bool isInterleaved = false)
This method will be called once FSR has won the competition with MPI 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.
Normally the optional `isInterleaved` argument would not be of interest. It can be used to separate resonance decays, false, from the interleaved evolution together with MPI and ISR, true. More precisely, it separates calls to the `timesDecPtr` and the `timesPtr` instances.

virtual bool TimeShower::rescatterPropogateRecoil( Event& event, Vec4& pNew)
This method is only called if rescattering is switched on in the description of multiparton interactions. It then propagates a recoil from a timelike branching to internal lines that connect systems. As for `rescatterUpdate` above, this is not likely to be of interest to most implementors of new showers.

int TimeShower::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 occurred, 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`).

virtual void TimeShower::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 virtual methods in the following. 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.

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

virtual SpaceShower::~SpaceShower()
Also the destructor does not need to do anything.

void SpaceShower::initPtr(Info* infoPtrIn, Settings* settingsPtrIn, ParticleData* particleDataPtrIn, Rndm* rndmPtrIn, CoupSM* coupSMPtrIn, PartonSystems* partonSystemsPtrIn, UserHooks* userHooksPtrIn, MergingHooks* mergingHooksPtrIn = 0)
This method only imports pointers to standard facilities, and is not virtual.

virtual void SpaceShower::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.

virtual bool SpaceShower::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 double-count. The `SpaceShower:pTmaxMatch` switch allows you to force the behaviour, or else to program your own logic. The current default implementation limits pT whenever the final state contains a quark (except top), gluon or photon, since then the danger of double-counting 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.

virtual double SpaceShower::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.

virtual void SpaceShower::prepare( int iSys, Event& event, bool limitPTmaxIn = true, double pTfirstTrialIn = 1e9)
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 MPI, since there an unlimited pT for sure would lead to double-counting.
The `pTfirstTrialIn` parameter is tentatively used in the new weak-shower part of the code. It is related to the option to shift the weak-shower evolution scale up by (some factor times) the gauge boson mass, without this backfiring on the QCD/QED evolution. Still under evolution, and normally not used.

virtual void SpaceShower::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.

virtual double SpaceShower::pTnext( Event& event, double pTbegAll, double pTendAll, int nRadIn = -1, bool isFirstTrial = false)
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 MPI 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 MPI in the interleaved evolution, but this is an implementation 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 MPI, ISR and FSR still have to be compared to assign the winner.
The 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.
The `isFirstTrialIn` is `true` only for the very first emission, and can then optionally be used to switch the weak-shower evolution scale up by (some factor times) the gauge boson mass.

virtual bool SpaceShower::branch( Event& event)
This method will be called once ISR has won the competition with MPI and FSR 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. Also a complete restart of the parton-level description may be necessary, see `doRestart()` below.

int SpaceShower::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 occurred, 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`).

bool SpaceShower::doRestart()
This method is not virtual. If `branch(...)` above fails to construct a branching, and the conditions are such that the whole parton-level description should be restarted, then it should return true, else not. Currently only the rescattering description can give this kind of failures, and so the internal `rescatterFail` boolean must be set true when this should happen, and else false.

virtual void SpaceShower::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.