1 <chapter name="Implement New Showers">
3 <h2>Implement New Showers</h2>
5 In case you want to replace the PYTHIA initial- and final-state
6 showers by your own, it is possible but not trivial. The point is
7 that multiple interactions (MI), initial-state radiation (ISR) and
8 final-state radiation (FSR) in general appear in one single
9 interleaved sequence of decreasing <ei>pT</ei> values. Therefore
10 shower replacements would have to be able to play the game by such
11 rules, as we will outline further below. Of course, this still
12 leaves the field open exactly how to define what to mean by
13 <ei>pT</ei>, how to handle recoil effects, how the colour flow is
14 affected, and so on, so there is certainly room for alternative
15 showers. A first example of a shower implemented within the PYTHIA
16 context is <a href="http://home.fnal.gov/~skands/vincia/">VINCIA</a>,
17 which however so far only handles FSR.
20 For the moment we assume you want to keep the MI part of the story
21 unchanged, and make use of the existing beam-remnants (BR) machinery.
22 If you want to replace both MI, ISR, FSR and BR then you had better
23 replace the whole <code>PartonLevel</code> module of the code.
24 If, in addition, you want to produce your own hard processes,
25 then you only need the
26 <aloc href="HadronLevelStandalone">hadron-level standalone</aloc>
27 part of the machinery.
30 In order to write replacement codes for ISR and/or FSR it is useful
31 to be aware of which information has to be shared between the
32 different components, and which input/output structure is required
33 of the relevant methods. For details, nothing beats studying the
34 existing code. However, here we provide an overview, that should
35 serve as a useful introduction.
38 It should be noted that we here primarily address the problem in
39 its full generality, with interleaved MI, ISR and FSR. There exists
40 an option <code>TimeShower:interleave = off</code> where only
41 MI and ISR would be interleaved and FSR be considered after these
42 two, but still before BR. Most of the aspects described here would
43 apply also for that case. By contrast, resonance decays are only
44 considered after all the four above components, and timelike
45 showers in those decays would never be interleaved with anything
46 else, so are much simpler to administrate.
49 Therefore the <code><aloc href="ProgramFlow">
50 pythia.setShowerPtr( timesDecPtr, timesPtr, spacePtr)</aloc></code>
51 method allows two separate pointers to be set to instances of
52 derived <code>TimeShower</code> classes. The first is only required
53 to handle decays, say of <ei>Z^0</ei> or <ei>Upsilon</ei>, with no
54 dependence on beam remnants or ISR. The second, as well as
55 <code>spacePtr</code>, has to handle the interleaved evolution of MI,
56 ISR and FSR. Therefore you are free to implement only the first, and
57 let the PYTHIA default showers take care of the latter two. But, if
58 you wanted to, you could also set <code>timesDecPtr = 0</code> and
59 only provide a <code>timesPtr</code>, or only a <code>spacePtr</code>.
60 If your timelike shower does both cases, the first two pointers
61 can agree. The only tiny point to take into account then is that
62 <code>init( beamAPtr, beamBPtr)</code> is called twice, a first time
63 to <code>timesDecPtr</code> with beam pointers 0, and a second time
64 to <code>timesPtr</code> with nonvanishing beam pointers.
66 <h3>The event record and associated information</h3>
68 Obviously the main place for sharing information is the event
69 record, specifically the <code>Event event</code> member of
70 <code>Pythia</code>, passed around as a reference. It is
71 assumed you already studied how it works, so here we only draw
72 attention to a few aspects of special relevance.
75 One basic principle is that existing partons should not be
76 overwritten. Instead new partons should be created, even when a
77 parton only receives a slightly shifted momentum and for the rest
78 stays the same. Such "carbon copies" by final-state branchings
79 should be denoted by both daughter indices of the original parton
80 pointing to the copy, and both mother indices of the copy to the
81 original. If the copy instead is intended to represent an earlier
82 step, e.g. in ISR backwards evolution, the role of mothers and
83 daughters is interchanged. The
84 <code>event.copy( iCopy, newStatus)</code>
85 routine can take care of this tedious task; the sign of
86 <code>newStatus</code> tells the program which case to assume.
89 To make the event record legible it is essential that the
90 <aloc href="ParticleProperties">status codes</aloc>
91 are selected appropriately to represent the reason why each new
92 parton is added to the record. Also remember to change the
93 status of a parton to be negative whenever an existing parton
94 is replaced by a set of new daughter partons.
97 Another important parton property is <code>scale()</code>,
98 which does not appear in the normal event listing, but only
99 if you use the extended
100 <code>Event:listScaleAndVertex = on</code> option.
101 This property is supposed to represent the production scale
102 (in GeV) of a parton. In the current FSR and ISR algorithms
103 it is used to restrict from above the allowed <ei>pT</ei>
104 values for branchings of this particular parton.
105 Beam remnants and other partons that should not radiate are
109 Auxiliary to the event record proper is the
110 <code><aloc href="AdvancedUsage">PartonSystems</aloc></code>
111 class, that keep track of which partons belong together in the
112 same scattering subsystem. This information must be kept up-to-date
113 during the shower evolution.
116 For initial-state showers it is also necessary to keep track of
117 the partonic content extracted from the beams. This information
119 <code><aloc href="AdvancedUsage">BeamParticle</aloc></code>
122 <h3>The TimeShower interface</h3>
124 If you want to replace the <code>TimeShower</code> class this would
125 involve replacing the virtual methods among the following ones.
127 <method name="TimeShower::TimeShower()">
128 The constructor does not need to do anything.
131 <method name="virtual TimeShower::~TimeShower()">
132 The destructor does not need to do anything.
135 <method name="void TimeShower::initPtr(Info* infoPtr, Settings* settingsPtr,
136 ParticleData* particleDataPtr, Rndm* rndmPtr, CoupSM* coupSMPtr,
137 PartonSystems* partonSystemsPtr, UserHooks* userHooksPtr)">
138 This method only imports pointers to standard facilities,
142 <method name="virtual void TimeShower::init( BeamParticle* beamAPtrIn = 0,
143 BeamParticle* beamBPtrIn = 0)">
144 You have to store your local copy of the pointers to these objects,
145 since they have to be used during the generation, as explained above.
146 The pointers could be zero; e.g. a local copy of <code>TimeShower</code>
147 is created to handle showers in decays such as <ei>Upsilon -> q qbar</ei>
148 from inside the <code>ParticleDecays class</code>. This is also the
149 place to do initialization of whatever parameters you plan to use,
150 e.g. by reading in them from a user-accessible database like the
151 <code>Settings</code> one.
154 <method name="virtual bool TimeShower::limitPTmax( Event& event,
155 double Q2Fac = 0., double Q2Ren = 0.)">
156 The question is whether the FSR should be allowed to occur at larger
157 scales than the hard process it surrounds. This is process-dependent,
158 as illustrated below for the the analogous
159 <code>SpaeShower::limitPTmax(...)</code> method, although the two
160 kinds of radiation need not have to be modelled identically.
161 The <code>TimeShower:pTmaxMatch</code> switch allows you to force the
162 behaviour among three options, but you may replace by your own logic.
163 <br/>The internal PYTHIA implementation also allows intermediate options,
164 where emissions can go up to the kinematical limit but be dampened above
165 the factorization or renormalization scale. Therefore the (square of the)
166 latter two are provided as optional input parameters.
169 <method name="double TimeShower::enhancePTmax()">
170 Relative to the default <ei>pT_max</ei> evolution scale of the process,
171 it may still be convenient to vary the matching slightly for the hardest
172 interaction in an event, to probe the sensitivity to such details. The
173 base-class implementation returns the value of the
174 <code>TimeShower:pTmaxFudge</code> parameter.
177 <method name="virtual int TimeShower::shower( int iBeg, int iEnd,
178 Event& event, double pTmax, int nBranchMax = 0)">
179 This is an all-in-one call for shower evolution, and as such cannot be
180 used for the normal interleaved evolution, where only the routines below
181 are used. It also cannot be used in resonance decays that form part of
182 the hard process, since there the
183 <aloc href="UserHooks">user hooks</aloc> insert a potential
184 veto step. Currently this routine is therefore only used in the
185 hadron-level decays, e.g. <ei>Upsilon -> g g g</ei>.
186 <br/><code>iBeg</code> and <code>iEnd</code> is the position of the
187 first and last parton of a separate system, typically produced by a
188 resonance decay. Such a system only evolves in isolation, and in
189 particular does not relate to the beams.
190 <br/>The <code>pTmax</code> value sets the maximum scale for evolution,
191 but normally you would restrict that further for each individual
192 parton based on its respective scale value.
193 <br/>The <code>nBranchMax</code> value, if positive, gives the maximum
194 number of allowed branchings in the call, as useful for matching studies.
195 <br/>The routine is expected to return the number of FSR branchings that
196 were generated, but only for non-critical statistics purposes.
197 <br/>Since the real action typically is delegated to the routines
198 below, it may well be that the existing code need not be replaced.
201 <method name="double TimeShower::pTLastInShower()">
202 Can be used to return the <ei>pT</ei> evolution scale of the last
203 branching in the cascade generated with the above
204 <code>shower(...)</code> method. Is to be set in the internal
205 <code>pTLastInShower</code> variable, and should be 0 if there
206 were no branchings. Can be useful for matching studies.
209 <method name="virtual void TimeShower::prepare( int iSys, Event& event)">
210 This method is called immediately after a new interaction (or the
211 products of a resonance decay) has been added, and should then be used
212 to prepare the subsystem of partons for subsequent evolution. In
213 the current code this involves identifying all colour and charge
214 dipole ends: the position of radiating and recoiling partons, maximum
215 <ei>pT</ei> scales, possible higher-order matrix elements matchings
217 <br/>The <code>iSys</code> parameter specifies which parton system
218 is to be prepared. It is used to extract the set of partons to be
219 treated, with rules as described in the above section on subsystems.
220 Specifically, the first two partons represent the incoming state,
221 or are 0 for resonance decays unrelated to the beams, while the
222 rest are not required to be in any particular order.
225 <method name="virtual void TimeShower::rescatterUpdate( int iSys,
227 This method is called immediately after rescattering in the description
228 of multiple interactions. Thus the information on one or several
229 systems is out-of-date, while that of the others is unchanged.
230 We do not provide the details here, since we presume few implementors
231 of new showers will want to touch the technicalities involved
232 in obtaining a description of rescattering.
235 <method name="virtual void TimeShower::update( int iSys, Event& event)">
236 This method is called immediately after a spacelike branching in the
237 <code>iSys</code>'th subsystem. Thus the information for that system is
238 out-of-date, while that of the others is unchanged. If you want, you are
239 free to throw away all information for the affected subsystem and call
240 <code>prepare( iSys, event)</code> to create new one. Alternatively
241 you may choose only to update the information that has changed.
244 <method name="virtual double TimeShower::pTnext( Event& event,
245 double pTbegAll, double pTendAll)">
246 This is the main driver routine for the downwards evolution. A new
247 <ei>pT</ei> is to be selected based on the current information set up
248 by the routines above, and along with that a branching parton or dipole.
249 The <code>pTbegAll</code> scale is the maximum scale allowed, from which
250 the downwards evolution should be begun (usually respecting the maximum
251 scale of each individual parton). If no emission is found above
252 <code>pTendAll</code> (and above the respective shower cutoff scales)
253 then <code>0.</code> should be returned and no emissions will be allowed.
254 Both scales can vary from one event to the next: if a scale has
255 already been selected for MI or ISR it makes no sense to look for
256 a scale smaller than that from FSR, since it would not be able to
257 compete, so <code>pTendAll</code> is set correspondingly. As it happens,
258 FSR is tried before ISR and MI in the interleaved evolution,
259 but this is an implementational detail that could well change.
260 <br/>Typically the implementation of this routine would be to set
261 up a loop over all possible radiating objects (dipoles, dipole ends, ...),
262 for each pick its possible branching scale and then pick the one
263 with largest scale as possible winner. At this stage no branching should
264 actually be carried out, since MI, ISR and FSR still have to be compared
265 to assign the winner.
268 <method name="virtual bool TimeShower::branch( Event& event,
269 bool isInterleaved = false)">
270 This method will be called once FSR has won the competition with
271 MI and ISR to do the next branching. The candidate branching found
272 in the previous step should here be carried out in full. The
273 pre-branching partons should get a negative status code and new
274 replacement ones added to the end of the event record. Also the
275 subsystem information should be updated, and possibly also the
277 <br/>Should some problem be encountered in this procedure, e.g. if
278 some not-previously-considered kinematics requirement fails, it is
279 allowed to return <code>false</code> to indicate that no branching
280 could be carried out.
281 <br/>Normally the optional <code>isInterleaved</code> argument would
282 not be of interest. It can be used to separate resonance decays, false,
283 from the interleaved evolution together with MI and ISR, true.
284 More precisely, it separates calls to the <code>timesDecPtr</code>
285 and the <code>timesPtr</code> instances.
288 <method name="virtual bool TimeShower::rescatterPropogateRecoil(
289 Event& event, Vec4& pNew)">
290 This method is only called if rescattering is switched on in the
291 description of multiple interactions. It then propagates a recoil
292 from a timelike branching to internal lines that connect systems.
293 As for <code>rescatterUpdate</code> above, this is not likely to be
294 of interest to most implementors of new showers.
297 <method name="int TimeShower::system()">
298 This method is not virtual. If a branching is constructed by the
299 previous routine this tiny method should be able to return the number
300 of the selected subsystem <code>iSysSel</code> where it occured,
301 so that the spacelike shower can be told which system to update,
302 if necessary. Therefore <code>iSysSel</code> must be set in
303 <code>branch</code> (or already in <code>pTnext</code>).
306 <method name="virtual void TimeShower::list( ostream& os = cout)">
307 This method is not at all required. In the current implementation it
308 outputs a list of all the dipole ends, with information on the
309 respective dipole. The routine is not called anywhere in the public
310 code, but has been inserted at various places during the
311 development/debug phase.
314 <h3>The SpaceShower interface</h3>
316 If you want to replace the <code>SpaceShower</code> class this would
317 involve replacing the virtual methods in the following. You will find
318 that much of the story reminds of <code>TimeShower</code> above, and
319 actually some cut-and-paste of text is involved. In some respects the
320 description is simpler, since there are no special cases for resonance
321 decays and non-interleaved evolution. Thus there is no correspondence
322 to the <code>TimeShower::shower(...)</code> routine.
324 <method name="SpaceShower::SpaceShower()">
325 The constructor does not need to do anything.
328 <method name="virtual SpaceShower::~SpaceShower()">
329 Also the destructor does not need to do anything.
332 <method name="void SpaceShower::initPtr(Info* infoPtr,
333 Settings* settingsPtr, ParticleData* particleDataPtr, Rndm* rndmPtr,
334 PartonSystems* partonSystemsPtr, UserHooks* userHooksPtr)">
335 This method only imports pointers to standard facilities,
339 <method name="virtual void SpaceShower::init( BeamParticle* beamAPtrIn,
340 BeamParticle* beamBPtrIn)">
341 You have to store your local copy of the pointers to these objects,
342 since they have to be used during the generation, as explained above.
343 This is also the place to do initialization of whatever parameters
344 you plan to use, e.g. by reading in them from a user-accessible
345 database like the <code>Settings</code> one.
348 <method name="virtual bool SpaceShower::limitPTmax( Event& event,
349 double Q2Fac = 0., double Q2Ren = 0.)">
350 The question is whether the ISR should be allowed to occur at larger
351 scales than the hard process it surrounds. This is process-dependent.
352 For instance, if the hard process is <ei>Z^0</ei> production we know
353 that ISR should be allowed to go right up to the kinematical limit.
354 If it is a <ei>2 -> 2</ei> QCD process the ISR should not exceed
355 the scale of the hard process, since if so one would doublecount.
356 The <code>SpaceShower:pTmaxMatch</code> switch allows you to force the
357 behaviour, or else to program your own logic. The current default
358 implementation limits <ei>pT</ei> whenever the final state contains
359 a quark (except top), gluon or photon, since then the danger of
360 doublecounting is there. You may replace by your own logic, or leave as is.
361 <br/>The internal PYTHIA implementation also allows intermediate options,
362 where emissions can go up to the kinematical limit but be dampened above
363 the factorization or renormalization scale. Therefore the (square of the)
364 latter two are provided as optional input parameters.
367 <method name="virtual double SpaceShower::enhancePTmax()">
368 When the above method limits <ei>pT_max</ei> to the scale of the process,
369 it may still be convenient to vary the matching slightly for the hardest
370 interaction in an event, to probe the sensitivity to such details. The
371 base-class implementation returns the value of the
372 <code>SpaceShower:pTmaxFudge</code> parameter.
375 <method name="virtual void SpaceShower::prepare( int iSys, Event& event,
376 bool limitPTmax = true)">
377 This method is called immediately after a new interaction has been
378 added, and should then be used to prepare the subsystem of partons
379 for subsequent evolution. In the current code this involves identifying
380 the colour and charge dipole ends: the position of radiating and recoiling
381 partons, maximum <ei>pT</ei> scales, and possible higher-order matrix
382 elements matchings to apply. Depending on what you have in mind you may
383 choose to store slightly different quantities. You have to use the
384 subsystem information described above to find the positions of the two
385 incoming partons (and the outgoing ones) of the system, and from there
386 the scales at which they were produced.
387 <br/> The <code>limitPTmax</code> input agrees with the output of the
388 previous method for the hardest process, and is always true for
389 subsequent MI, since there an unlimited <ei>pT</ei> for sure
390 would lead to doublecounting.
393 <method name="virtual void SpaceShower::update( int iSys, Event& event)">
394 This method is called immediately after a timelike branching in the
395 <code>iSys</code>'th subsystem. Thus the information for that system may
396 be out-of-date, and to be updated. For the standard PYTHIA showers
397 this routine does not need to do anything, but that may be different
398 in another implementation.
401 <method name="virtual double SpaceShower::pTnext( Event& event,
402 double pTbegAll, double pTendAll, int nRadIn = -1)">
403 This is the main driver routine for the downwards evolution. A new
404 <ei>pT</ei> is to be selected based on the current information set up
405 by the routines above, and along with that a branching parton or dipole.
406 The <code>pTbegAll</code> scale is the maximum scale allowed, from which
407 the downwards evolution should be begun (usually respecting the maximum
408 scale of each individual parton). If no emission is found above
409 <code>pTendAll</code> (and above the respective shower cutoff scales)
410 then <code>0.</code> should be returned and no emissions will be allowed.
411 Both scales can vary from one event to the next: if a scale has
412 already been selected for MI or ISR it makes no sense to look for
413 a scale smaller than that from FSR, since it would not be able to
414 compete, so <code>pTendAll</code> is set correspondingly. As it happens,
415 FSR is tried before ISR and MI in the interleaved evolution,
416 but this is an implementational detail that could well change.
417 <br/>Typically the implementation of this routine would be to set
418 up a loop over all possible radiating objects (dipoles, dipole ends, ...),
419 for each pick its possible branching scale and then pick the one
420 with largest scale as possible winner. At this stage no branching should
421 actually be carried out, since MI, ISR and FSR still have to be compared
422 to assign the winner.
423 <br/>The final input <code>nRadIn</code> provides the total number of
424 ISR and FSR emissions already generated in the event, and so allows a
425 special treatment for the very first emission, if desired.
428 <method name="virtual bool SpaceShower::branch( Event& event)">
429 This method will be called once FSR has won the competition with
430 MI and ISR to do the next branching. The candidate branching found
431 in the previous step should here be carried out in full. The
432 pre-branching partons should get a negative status code and new
433 replacement ones added to the end of the event record. Also the
434 subsystem information should be updated, and possibly also the
436 <br/>Should some problem be encountered in this procedure, e.g. if
437 some not-previously-considered kinematics requirement fails, it is
438 allowed to return <code>false</code> to indicate that no branching
439 could be carried out. Also a complete restart of the parton-level
440 description may be necessary, see <code>doRestart()</code> below.
443 <method name="int SpaceShower::system()">
444 This method is not virtual. If a branching is constructed by the
445 previous routine this tiny method should be able to return the number
446 of the selected subsystem <code>iSysSel</code> where it occured,
447 so that the spacelike shower can be told which system to update,
448 if necessary. Therefore <code>iSysSel</code> must be set in
449 <code>branch</code> (or already in <code>pTnext</code>).
452 <method name="bool SpaceShower::doRestart()">
453 This method is not virtual. If <code>branch(...)</code> above fails
454 to construct a branching, and the conditions are such that the whole
455 parton-level description should be restarted, then it should return
456 true, else not. Currently only the rescattering description can give
457 thi kind of failures, and so the internal <code>rescatterFail</code>
458 boolean must be set true when this should happen, and else false.
461 <method name="virtual void SpaceShower::list( ostream& os = cout)">
462 This method is not at all required. In the current implementation it
463 outputs a list of all the dipole ends, with information on the
464 respective dipole. The routine is not called anywhere in the public
465 code, but has been inserted at various places during the
466 development/debug phase.
471 <!-- Copyright (C) 2010 Torbjorn Sjostrand -->