- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / xmldoc / ImplementNewShowers.xml
1 <chapter name="Implement New Showers">
2
3 <h2>Implement New Showers</h2>
4
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. 
18
19 <p/>
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. 
28
29 <p/>
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.
36
37 <p/>
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.
47
48 <p/>
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.  
65
66 <h3>The event record and associated information</h3>  
67
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.
73
74 <p/>
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.
87
88 <p/>
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.
95
96 <p/>
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 
106 assigned scale 0.   
107
108 <p/>
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.
114
115 <p/>
116 For initial-state showers it is also necessary to keep track of
117 the partonic content extracted from the beams. This information
118 is stored in the 
119 <code><aloc href="AdvancedUsage">BeamParticle</aloc></code>  
120 class.  
121
122 <h3>The TimeShower interface</h3>
123
124 If you want to replace the <code>TimeShower</code> class this would
125 involve replacing the virtual methods among the following ones.
126
127 <method name="TimeShower::TimeShower()">
128 The constructor does not need to do anything.
129 </method>
130
131 <method name="virtual TimeShower::~TimeShower()">
132 The destructor does not need to do anything.
133 </method>
134
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, 
139 and is not virtual.
140 </method>
141
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.
152 </method>
153
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.  
167 </method>
168
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.
175 </method>
176
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.
199 </method>
200
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. 
207 </method>
208
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 
216 to apply, and so on. 
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.
223 </method>
224
225 <method name="virtual void TimeShower::rescatterUpdate( int iSys, 
226 Event& event)">
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.
233 </method>
234
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.
242 </method>
243
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.
266 </method>
267
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
276 beams. 
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.
286 </method>
287
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.
295 </method>
296
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>).  
304 </method>
305
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.
312 </method>
313    
314 <h3>The SpaceShower interface</h3>
315
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.  
323
324 <method name="SpaceShower::SpaceShower()">
325 The constructor does not need to do anything.
326 </method>
327
328 <method name="virtual SpaceShower::~SpaceShower()">
329 Also the destructor does not need to do anything.
330 </method>
331
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, 
336 and is not virtual.
337 </method>
338
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.
346 </method>
347
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.  
365 </method>
366
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.
373 </method>
374
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.
391 </method>
392
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.
399 </method>
400
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.
426 </method>
427
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
435 beams. 
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.  
441 </method>
442
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>).  
450 </method>
451
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.
459 </method>
460
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.
467 </method>
468    
469 </chapter>
470
471 <!-- Copyright (C) 2010 Torbjorn Sjostrand -->