using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / xmldoc / ImplementNewShowers.xml
CommitLineData
5ad4eb21 1<chapter name="Implement New Showers">
2
3<h2>Implement New Showers</h2>
4
5In case you want to replace the PYTHIA initial- and final-state
6showers by your own, it is possible but not trivial. The point is
7that multiple interactions (MI), initial-state radiation (ISR) and
8final-state radiation (FSR) in general appear in one single
9interleaved sequence of decreasing <ei>pT</ei> values. Therefore
10shower replacements would have to be able to play the game by such
11rules, as we will outline further below. Of course, this still
12leaves 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
14affected, and so on, so there is certainly room for alternative
15showers.
16
17<p/>
18For the moment we assume you want to keep the MI part of the story
19unchanged, and make use of the existing beam-remnants (BR) machinery.
20If you want to replace both MI, ISR, FSR and BR then you had better
21replace the whole <code>PartonLevel</code> module of the code.
22If, in addition, you want to produce your own hard processes,
23then you only need the
24<aloc href="HadronLevelStandalone">hadron-level standalone</aloc>
25part of the machinery.
26
27<p/>
28In order to write replacement codes for ISR and/or FSR it is useful
29to be aware of which information has to be shared between the
30different components, and which input/output structure is required
31of the relevant methods. For details, nothing beats studying the
32existing code. However, here we provide an overview, that should
33serve as a useful introduction.
34
35<p/>
36It should be noted that we here primarily address the problem in
37its full generality, with interleaved MI, ISR and FSR. There exists
38an option <code>TimeShower:interleave = off</code> where only
39MI and ISR would be interleaved and FSR be considered after these
40two, but still before BR. Most of the aspects described here would
41apply also for that case. By contrast, resonance decays are only
42considered after all the four above components, and timelike
43showers in those decays would never be interleaved with anything
44else, so are much simpler to administrate.
45
46<p/>
47Therefore the <aloc href="ProgramFlow">
48<code>pythia.setShowerPtr( timesDecPtr, timesPtr, spacePtr)</code></aloc>
49method allows two separate pointers to be set to instances of
50derived <code>TimeShower</code> classes. The first is only required
51to handle decays, say of <ei>Z^0</ei> or <ei>Upsilon</ei>, with no
52dependence on beam remnants or ISR. The second, as well as
53<code>spacePtr</code>, has to handle the interleaved evolution of MI,
54ISR and FSR. Therefore you are free to implement only the first, and
55let the PYTHIA default showers take care of the latter two. But, if
56you wanted to, you could also set <code>timesDecPtr = 0</code> and
57only provide a <code>timesPtr</code>, or only a <code>spacePtr</code>.
58If your timelike shower does both cases, the first two pointers
59can agree. The only tiny point to take into account then is that
60<code>init( beamAPtr, beamBPtr)</code> is called twice, a first time
61to <code>timesDecPtr</code> with beam pointers 0, and a second time
62to <code>timesPtr</code> with nonvanishing beam pointers.
63
64<h3>The event record</h3>
65
66Obviously the main place for sharing information is the event
67record, specifically the <code>Event event</code> member of
68<code>Pythia</code>, passed around as a reference. It is
69assumed you already studied how it works, so here we only draw
70attention to a few aspects of special relevance.
71
72<p/>
73One basic principle is that existing partons should not be
74overwritten. Instead new partons should be created, even when a
75parton only receives a slightly shifted momentum and for the rest
76stays the same. Such "carbon copies" by final-state branchings
77should be denoted by both daughter indices of the original parton
78pointing to the copy, and both mother indices of the copy to the
79original. If the copy instead is intended to represent an earlier
80step, e.g. in ISR backwards evolution, the role of mothers and
81daughters is interchanged. The
82<code>event.copy( iCopy, newStatus)</code>
83routine can take care of this tedious task; the sign of
84<code>newStatus</code> tells the program which case to assume.
85
86<p/>
87To make the event record legible it is essential that the
88<aloc href="ParticleProperties">status codes</aloc>
89are selected appropriately to represent the reason why each new
90parton is added to the record. Also remember to change the
91status of a parton to be negative whenever an existing parton
92is replaced by a set of new daughter partons.
93
94<p/>
95Another important parton property is <code>scale()</code>,
96which does not appear in the normal event listing, but only
97if you use the extended
98<code>Event:listScaleAndVertex = on</code> option.
99This property is supposed to represent the production scale
100(in GeV) of a parton. In the current FSR and ISR algorithms
101it is used to restrict from above the allowed <ei>pT</ei>
102values for branchings of this particular parton.
103Beam remnants and other partons that should not radiate are
104assigned scale 0.
105
106<h3>The subsystems</h3>
107
108One aspect that complicates administration is that an event
109can contain several subsystems, consisting of one MI and its
110associated ISR and FSR, which to first approximation are assumed
111to evolve independently, but to second are connected by the
112interleaved evolution. The partons of a given subsystem
113therefore do not have to be stored consecutively. Some simple
114vectors put inside the event record can be used to keep track of
115the current position, i.e. index <code>iPos</code> for parton
116<code>event[iPos]</code>, of all partons of all systems.
117
118<p/>
119The number of systems is given by <code>sizeSystems()</code>,
120with systems numbered 0 through <code>sizeSystems() - 1</code>.
121The size of a given subsystem <code>iSys</code> is given by
122<code>sizeSystem(iSys)</code>, with analogous numbering.
123The method <code>getInSystem( iSys, iMem)</code> returns the
124position <code>iPos</code> of the <code>iMem</code>'th member of
125the <code>iSys</code>'th system. For each system, the slots
126<code>iMem =</code> 0 and 1 are intended for the incoming partons
127from beam A and B, respectively. If there are no beams, such as
128in resonance decays, these should be assigned value 0.
129Slots 2 onwards are intended for all the outgoing partons. These
130latter partons are not guaranteed to appear in any particular order.
131
132<p/>
133New systems are created from the hard process and by the MI,
134not from any of the other components, by the <code>newSystem()</code>
135method, which returns the <code>iSys</code> code of the new system.
136Both FSR and ISR modify the position of partons, however.
137Since an FSR or ISR branching typically implies a new state with
138one more parton than before, the method
139<code>addToSystem( iSys, iPos)</code> appends a new slot to the
140given system and stores the <code>iPos</code> value there. Furthermore,
141in a branching, several existing partons may also be moved to new
142slots. The method <code>setInSystem( iSys, iMem, iPos)</code>
143replaces the current <code>iPos</code> value in the given slot
144by the provided one. If you do not know the <code>iMem</code> value,
145<code>replaceInSystem( iSys, iPosOld, iPosNew)</code>
146will replace any ocurence of <code>iPosOld</code> by
147<code>iPosNew</code> wherever it is found in the
148<code>iSys</code>'th system. In a FSR <ei>1 -> 2</ei> branching
149it is irrelevant which parton position you let overwrite the
150existing slot and which is added to the end of the system.
151
152<p/>
153Finally, <code>clearSystems()</code> empties the information, and
154<code>listSystems()</code> provides a simple listing of all the
155<code>iPos</code> values by system, intended for debugging purposes only.
156
157<p/>
158The system information must be kept up-to-date. The MI component only
159writes, but both ISR, FSR and BR make extensive use of the existing
160information. As an example, the introduction of primordial <ei>kT</ei>
161in the beam remnants will fail if the information on which
162final-state partons belong to which system is out-of-date.
163
164<p/>
165Currently the system information is kept throughout the continued
166history of the event. Resonance decays create new systems, appended
167to the existing ones. This could be useful during the hadronization
168stage, to collect the partons that belong to a resonace with
169preserved mass when a small string collapses to one particle,
170but is not yet used for that.
171
172<h3>The beams</h3>
173
174The different subsystems are tied together by them sharing the same
175initial beam particles, and thereby being restricted by energy-momentum
176and flavour conservation issues. The information stored in the two
177beam particles, here called <code>beamA</code> and <code>beamB</code>,
178is therefore as crucial to keep correct as the above subsystem list.
179
180<p/>
181Both beam objects are of the <code>BeamParticle</code> class.
182Each such object contains a vector with the partons extracted from it.
183The number of such partons, <code>beamX.size()</code> (X = A or B),
184of course is the same as the above number of subsystems in the event
185record. (The two diverge at the BR step, where further beam remnants
186are added to the beams without corresponding to new subsystems.)
187The individual partons are accessed by an overloaded indexing
188operator to a vector of <code>ResolvedParton</code> objects. The
189<code>iPos()</code> property corresponds to the <code>iPos</code>
190one above, i.e. providing the position in the main event record of
191a parton. In particular,
192<code>beamA[iSys].iPos() = event.getInSystem( iSys, 0)</code> and
193<code>beamB[iSys].iPos() = event.getInSystem( iSys, 1)</code>.
194Whereas thus the indices of the two incoming partons to a subsystem
195are stored in two places, the ones of the outgoing partons only
196appear in the system part of the <code>Event</code> class.
197
198<p/>
199Just as the subsystems in <code>Event</code> must be updated, so must
200the information in the two <code>BeamParticle</code>'s, e.g. with methods
201<code>beamX[iSys].iPos( iPosIn)</code> when an incoming parton is
202replaced by a new one in line <code>iPosIn</code>. Furthermore the new
203parton identity should be set by <code>beamX[iSys].id( idIn)</code>
204and the new <ei>x</ei> energy-momentum fraction by
205<code>beamX[iSys].x( xIn)</code>. The three can be combined in one go
206by <code>beamX[iSys].update( iPosIn, idIn, xIn)</code>.
207
208<p/>
209To be specific, it is assumed that, at each step, the two incoming
210partons are moving along the <ei>+-z</ei> axis and are massless.
211Since the event is constructed in the c.m. frame of the incoming
212beams this implies that <ei>x = 2 E / E_cm</ei>.
213If the <ei>x</ei> values are not defined accordingly or not kept
214up-to-date the BR treatment will not conserve energy-momentum.
215
216<p/>
217In return, the <code>BeamParticle</code> objects give access to some
218useful methods. The <code>beamX.xf( id, x, Q2)</code> returns the
219standard PDF weight <ei>x f_id(x, Q^2)</ei>. More intererstingly,
220<code>beamX.xfISR( iSys, id, x, Q2)</code> returns the modified weight
221required when several subsystems have to share the energy and flavours.
222Thus <code>iSys</code> is added as an extra argument, and the momentum
223already assigned to the other subsystems is not available for evolution,
224i.e. the maximal <ei>x</ei> is correspondingly smaller than unity.
225Also flavour issues are handled in a similar spirit.
226
227<p/>
228An additional complication is that a parton can be either valence or
229sea, and in the latter case the BR treatment also distinguishes
230companion quarks, i.e. quark-antiquark pairs that are assumed to
231come from the same original <ei>g -> q qbar</ei> branching, whether
232perturbative or not. This can be queried either with the
233<code>beamX[iSys].companion()</code> method, for detailed information,
234or with the <code>beamX[iSys].isValence()</code>,
235<code>beamX[iSys].isUnmatched()</code> and
236<code>beamX[iSys].isCompanion()</code> metods for yes/no answers
237whether a parton is valence, unmatched sea or matched sea.
238This choice should affect the ISR evolution; e.g., a valence quark
239cannot be constructed back to come from a gluon.
240
241<p/>
242To keep this info up-to-date, the <code>beamX.pickValSeaComp()</code>
243method should be called whenever a parton of a new flavour has been
244picked in the ISR backwards evolution, but not if the flavour has not
245been changed, since then one should not be allowed to switch back and
246forth between the same quark being considered as valence or as sea.
247Since the <code>pickValSeaComp()</code> method makes use of the
248current parton-density values, it should be preceded by a call
249to <code>beamX.xfISR( iSys, id, x, Q2)</code>, where the values in
250the call are the now finally accepted ones for the newly-found mother.
251(Such a call is likely to have been made before, during the evolution,
252but is not likely to be the most recent one, i.e. still in memory, and
253therefore had better be redone.)
254
255<p/>
256Most of the issues in this section are related to the ISR algorithm,
257i.e. it is possible to write an FSR algorithm that is completely
258decoupled. Alternatively the FSR may change the position where an
259incoming parton is stored, or its assumed momentum, e.g. by recoil
260effects inside dipoles stretched from the scattered back to the
261incoming partons. In that case some of the methods above would have
262to be used also inside the FSR algorithm.
263
264<h3>The TimeShower interface</h3>
265
266If you want to replace the <code>TimeShower</code> class this would
267involve replacing the following methods.
268
269<method name="TimeShower()">
270The constructor does not need to do anything.
271</method>
272
273<method name="void init( BeamParticle* beamAPtrIn = 0,
274BeamParticle* beamBPtrIn = 0)">
275You have to store your local copy of the pointers to these objects,
276since they have to be used during the generation, as explained above.
277The pointers could be zero; e.g. a local copy of <code>TimeShower</code>
278is created to handle showers in decays such as <ei>Upsilon -> q qbar</ei>
279from inside the <code>ParticleDecays class</code>. This is also the
280place to do initialization of whatever parameters you plan to use,
281e.g. by reading in them from a user-accessible database like the
282<code>Settings</code> one.
283</method>
284
285<method name="double enhancePTmax()">
286Relative to the default <ei>pT_max</ei> evolution scale of the process,
287it may still be convenient to vary the matching slightly for the hardest
288interaction in an event, to probe the sensitivity to such details. The
289base-class implementation returns the value of the
290<code>TimeShower:pTmaxFudge</code> parameter.
291</method>
292
293<method name="int shower( int iBeg, int iEnd, Event& event, double pTmax)">
294This is an all-in-one call for shower evolution, and as such cannot be
295used for the normal interleaved evolution, where only the routines below
296are used. It also cannot be used in resonance decays that form part of
297the hard process, since there the
298<aloc href="UserHooks">user hooks</aloc> insert a potential
299veto step. Currently this routine is therefore only used in the
300hadron-level decays, e.g. <ei>Upsilon -> g g g</ei>.
301<br/><code>iBeg</code> and <code>iEnd</code> is the position of the
302first and last parton of a separate system, typically produced by a
303resonance decay. Such a system only evolves in isolation, and in
304particular does not relate to the beams.
305<br/>The <code>pTmax</code> value sets the maximum scale for evolution,
306but normally you would restrict that further for each individual
307parton based on its respective scale value.
308<br/>The routine is expected to return the number of FSR branchings that
309were generated, but only for non-critical statistics purposes.
310<br/>Since the real action typically is delegated to the routines
311below, it may well be that the existing code need not be replaced.
312</method>
313
314<method name="void prepare( int iSys, Event& event)">
315This method is called immediately after a new interaction (or the
316products of a resonance decay) has been added, and should then be used
317to prepare the subsystem of partons for subsequent evolution. In
318the current code this involves identifying all colour and charge
319dipole ends: the position of radiating and recoiling partons, maximum
320<ei>pT</ei> scales, possible higher-order matrix elements matchings
321to apply, and so on.
322<br/>The <code>iSys</code> parameter specifies which parton system
323is to be prepared. It is used to extract the set of partons to be
324treated, with rules as described in the above section on subsystems.
325Specifically, the first two partons represent the incoming state,
326or are 0 for resonance decays unrelated to the beams, while the
327rest are not required to be in any particular order.
328</method>
329
330<method name="void update( int iSys, Event& event)">
331This method is called immediately after a spacelike branching in the
332<code>iSys</code>'th subsystem. Thus the information for that system is
333out-of-date, while that of the others is unchanged. If you want, you are
334free to throw away all information for the affected subsystem and call
335<code>prepare( iSys, event)</code> to create new one. Alternatively
336you may choose only to update the information that has changed.
337</method>
338
339<method name="double pTnext( Event& event, double pTbegAll,
340double pTendAll)">
341This is the main driver routine for the downwards evolution. A new
342<ei>pT</ei> is to be selected based on the current information set up
343by the routines above, and along with that a branching parton or dipole.
344The <code>pTbegAll</code> scale is the maximum scale allowed, from which
345the downwards evolution should be begun (usually respecting the maximum
346scale of each individual parton). If no emission is found above
347<code>pTendAll</code> (and above the respective shower cutoff scales)
348then <code>0.</code> should be returned and no emissions will be allowed.
349Both scales can vary from one event to the next: if a scale has
350already been selected for MI or ISR it makes no sense to look for
351a scale smaller than that from FSR, since it would not be able to
352compete, so <code>pTendAll</code> is set correspondingly. As it happens,
353FSR is tried before ISR and MI in the interleaved evolution,
354but this is an implementational detail that could well change.
355<br/>Typically the implementation of this routine would be to set
356up a loop over all possible radiating objects (dipoles, dipole ends, ...),
357for each pick its possible branching scale and then pick the one
358with largest scale as possible winner. At this stage no branching should
359actually be carried out, since MI, ISR and FSR still have to be compared
360to assign the winner.
361</method>
362
363<method name="bool branch( Event& event)">
364This method will be called once FSR has won the competition with
365MI and ISR to do the next branching. The candidate branching found
366in the previous step should here be carried out in full. The
367pre-branching partons should get a negative status code and new
368replacement ones added to the end of the event record. Also the
369subsystem information should be updated, and possibly also the
370beams.
371<br/>Should some problem be encountered in this procedure, e.g. if
372some not-previously-considered kinematics requirement fails, it is
373allowed to return <code>false</code> to indicate that no branching
374could be carried out.
375</method>
376
377<method name="int system()">
378This method is not virtual. If a branching is constructed by the
379previous routine this tiny method should be able to return the number
380of the selected subsystem <code>iSysSel</code> where it occured,
381so that the spacelike shower can be told which system to update,
382if necessary. Therefore <code>iSysSel</code> must be set in
383<code>branch</code> (or already in <code>pTnext</code>).
384</method>
385
386<method name="void list( ostream& os = cout)">
387This method is not at all required. In the current implementation it
388outputs a list of all the dipole ends, with information on the
389respective dipole. The routine is not called anywhere in the public
390code, but has been inserted at various places during the
391development/debug phase.
392</method>
393
394<h3>The SpaceShower interface</h3>
395
396If you want to replace the <code>SpaceShower</code> class this would
397involve replacing the following methods. You will find that much of
398the story reminds of <code>TimeShower</code> above, and actually some
399cut-and-paste of text is involved. In some respects the description
400is simpler, since there are no special cases for resonance decays
401and non-interleaved evolution. Thus there is no correspondence to the
402<code>TimeShower::shower(...)</code> routine.
403
404<method name="SpaceShower()">
405The constructor does not need to do anything.
406</method>
407
408<method name="void init( BeamParticle* beamAPtrIn,
409BeamParticle* beamBPtrIn)">
410You have to store your local copy of the pointers to these objects,
411since they have to be used during the generation, as explained above.
412This is also the place to do initialization of whatever parameters
413you plan to use, e.g. by reading in them from a user-accessible
414database like the <code>Settings</code> one.
415</method>
416
417<method name="bool limitPTmax( Event& event, double Q2Fac = 0.,
418double Q2Ren = 0.)">
419The question is whether the ISR should be allowed to occur at larger
420scales than the hard process it surrounds. This is process-dependent.
421For instance, if the hard process is <ei>Z^0</ei> production we know
422that ISR should be allowed to go right up to the kinematical limit.
423If it is a <ei>2 -> 2</ei> QCD process the ISR should not exceed
424the scale of the hard process, since if so one would doublecount.
425The <code>SpaceShower:pTmaxMatch</code> switch allows you to force the
426behaviour, or else to program your own logic. The current implementation
427limits <ei>pT</ei> whenever the final state contains a quark (except top),
428gluon or photon, since then the danger of doublecounting is there.
429You may replace by your own logic, or leave as is.
430<br/>The internal PYTHIA implementation also allows intermediate options,
431where emissions can go up to the kinematical limit but be dampened above
432the factorization or renormalization scale. Therefore the (square of the)
433latter two are provided as optional input parameters.
434</method>
435
436<method name="double enhancePTmax()">
437When the above method limits <ei>pT_max</ei> to the scale of the process,
438it may still be convenient to vary the matching slightly for the hardest
439interaction in an event, to probe the sensitivity to such details. The
440base-class implementation returns the value of the
441<code>SpaceShower:pTmaxFudge</code> parameter.
442</method>
443
444<method name="void prepare( int iSys, Event& event,
445bool limitPTmax = true)">
446This method is called immediately after a new interaction has been
447added, and should then be used to prepare the subsystem of partons
448for subsequent evolution. In the current code this involves identifying
449the colour and charge dipole ends: the position of radiating and recoiling
450partons, maximum <ei>pT</ei> scales, and possible higher-order matrix
451elements matchings to apply. Depending on what you have in mind you may
452choose to store slightly different quantities. You have to use the
453subsystem information described above to find the positions of the two
454incoming partons (and the outgoing ones) of the system, and from there
455the scales at which they were produced.
456<br/> The <code>limitPTmax</code> input agrees with the output of the
457previous method for the hardest process, and is always true for
458subsequent MI, since there an unlimited <ei>pT</ei> for sure
459would lead to doublecounting.
460</method>
461
462<method name="void update( int iSys, Event& event)">
463This method is called immediately after a timelike branching in the
464<code>iSys</code>'th subsystem. Thus the information for that system may
465be out-of-date, and to be updated. For the standard PYTHIA showers
466this routine does not need to do anything, but that may be different
467in another implementation.
468</method>
469
470<method name="double pTnext(( Event& event, double pTbegAll,
471double pTendAll, int nRadIn = -1)">
472This is the main driver routine for the downwards evolution. A new
473<ei>pT</ei> is to be selected based on the current information set up
474by the routines above, and along with that a branching parton or dipole.
475The <code>pTbegAll</code> scale is the maximum scale allowed, from which
476the downwards evolution should be begun (usually respecting the maximum
477scale of each individual parton). If no emission is found above
478<code>pTendAll</code> (and above the respective shower cutoff scales)
479then <code>0.</code> should be returned and no emissions will be allowed.
480Both scales can vary from one event to the next: if a scale has
481already been selected for MI or ISR it makes no sense to look for
482a scale smaller than that from FSR, since it would not be able to
483compete, so <code>pTendAll</code> is set correspondingly. As it happens,
484FSR is tried before ISR and MI in the interleaved evolution,
485but this is an implementational detail that could well change.
486<br/>Typically the implementation of this routine would be to set
487up a loop over all possible radiating objects (dipoles, dipole ends, ...),
488for each pick its possible branching scale and then pick the one
489with largest scale as possible winner. At this stage no branching should
490actually be carried out, since MI, ISR and FSR still have to be compared
491to assign the winner.
492<br/>The final input <code>nRadIn</code> provides the total number of
493ISR and FSR emissions already generated in the event, and so allows a
494special treatment for the very first emission, if desired.
495</method>
496
497<method name="bool branch( Event& event)">
498This method will be called once FSR has won the competition with
499MI and ISR to do the next branching. The candidate branching found
500in the previous step should here be carried out in full. The
501pre-branching partons should get a negative status code and new
502replacement ones added to the end of the event record. Also the
503subsystem information should be updated, and possibly also the
504beams.
505<br/>Should some problem be encountered in this procedure, e.g. if
506some not-previously-considered kinematics requirement fails, it is
507allowed to return <code>false</code> to indicate that no branching
508could be carried out.
509</method>
510
511<method name="int system()">
512This method is not virtual. If a branching is constructed by the
513previous routine this tiny method should be able to return the number
514of the selected subsystem <code>iSysSel</code> where it occured,
515so that the spacelike shower can be told which system to update,
516if necessary. Therefore <code>iSysSel</code> must be set in
517<code>branch</code> (or already in <code>pTnext</code>).
518</method>
519
520<method name="void list( ostream& os = cout)">
521This method is not at all required. In the current implementation it
522outputs a list of all the dipole ends, with information on the
523respective dipole. The routine is not called anywhere in the public
524code, but has been inserted at various places during the
525development/debug phase.
526</method>
527
528</chapter>
529
530<!-- Copyright (C) 2008 Torbjorn Sjostrand -->