using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / xmldoc / SemiInternalProcesses.xml
1 <chapter name="Semi-Internal Processes">
2
3 <h2>Semi-Internal Processes</h2>
4
5 Normally users are expected to implement new processes via the
6 <aloc href="LesHouchesAccord">Les Houches Accord</aloc>. Then 
7 you do all flavour, colour and phase-space selection externally, 
8 before your process-level events are input for further processing
9 by PYTHIA. However, it is also possible to implement a 
10 new process in exactly the same way as the internal PYTHIA
11 ones, thus making use of the internal phase space selection machinery
12 to sample an externally provided cross-section expression. 
13 This page gives a brief summary how to do that. If you additionally 
14 want to introduce a new resonance species, with its own internal 
15 width calculations, you will find further instructions 
16 <aloc href="SemiInternalResonances">here</aloc>.
17
18 <p/>
19 Should you actually go ahead, it is strongly recommended to shop around 
20 for a similar process that has already been implemented, and to use that 
21 existing code as a template. Look for processes with the same combinations
22 of incoming flavours and colour flows, rather than the shape of the 
23 cross section itself. With a reasonable such match the task should be 
24 of medium difficulty, without it more demanding. 
25
26 <p/>
27 PYTHIA is rather good at handling the phase space of 
28 <ei>2 -> 1</ei> and <ei>2 -> 2</ei> processes, is more primitive for 
29 <ei>2 -> 3</ei> ones and does not at all address higher multiplicities. 
30 This limits the set of processes that you can implement in this 
31 framework. The produced particles may be resonances, however, so it is 
32 possible to end up with bigger "final" multiplicities through sequential 
33 decays, and to include further matrix-element weighting in those decays.   
34
35 <p/>
36 There are two steps involved in implementing a process:
37 <br/>1) writing a new class,  where the matrix elements are implemented, 
38 including information on incoming and outgoing flavours and colours, and 
39 <br/>2) making the process available.
40 <br/>We consider these two aspects in turn. An example where it all comes
41 together is found in <code>main25.cc</code>.
42
43 <h3>The Cross Section Class</h3> 
44
45 The matrix-element information has to be encoded in a new class.
46 The relevant code could either be put before the main program in the
47 same file, or be stored separately, e.g. in a matched pair
48 of <code>.h</code> and <code>.cc</code> files. The latter may be more
49 convenient, in particular if the cross sections are lengthy, or if you 
50 intend to build up your own little process library, but of course
51 requires that these additional files are correctly compiled and linked.
52
53 <p/>
54 The class has to be derived either from  
55 <code>Sigma1Process</code>, for <ei>2 -> 1</ei> processes, from 
56 <code>Sigma2Process</code>, for <ei>2 -> 2</ei> ones, or from 
57 <code>Sigma3Process</code>, for <ei>2 -> 3</ei> ones. (The 
58 <code>Sigma0Process</code> class is used for elastic, diffractive
59 and minimum-bias events, and is not recommended for use beyond that.) 
60 These are in their turn derived from the <code>SigmaProcess</code> 
61 base class. 
62
63 <p/>
64 The class can implement a number of methods. Some of these are 
65 compulsory, others strongly recommended, and the rest are to be
66 used only when the need arises to override the default behaviour. 
67 The methods are:
68
69 <p/>
70 A <b>constructor</b> for the derived class obviously must be available.
71 Here you are quite free to allow a list of arguments, to set
72 the parameters of your model, or even to create a set of closely
73 related but distinct processes. For instance, <ei>g g -> Q Qbar</ei>, 
74 <ei>Q = c</ei> or <ei>b</ei>, is only coded once, and then the 
75 constructor takes the quark code (4 or 5)  as argument, 
76 to allow the proper amount of differentiation. 
77
78 <p/>
79 A <b>destructor</b> is only needed if you plan to delete the process 
80 before the natural end of the run, and require some special behaviour 
81 at that point. If you call such a destructor you will leave a pointer 
82 dangling inside the <code>Pythia</code> object you gave it in to,
83 if that still exists. 
84
85 <method name="void initProc()"> 
86 is called once during initalization, and can then be used to set up
87 parameters, such as masses and couplings, and perform calculations 
88 that need not be repeated for each new event, thereby saving time. 
89 This method needs not be implemented, since in principle all 
90 calculations can be done in <code>sigmaHat</code> below.
91
92 <method name="void sigmaKin()">
93 is called once a kinematical configuration has been determined, but 
94 before the two incoming flavours are known. This routine can therefore 
95 be used to perform calculations that otherwise might have to be repeated 
96 over and over again in <code>sigmaHat</code> below. For instance 
97 a flavour-independent cross section calculation for a <ei>q g</ei> 
98 initial state would be repeated 20 times in <code>sigmaHat</code>, 
99 five times for the five quark flavours allowed in the incoming beams, 
100 times twice to include antiquarks, times twice since the (anti)quark 
101 could be in either of the two beams. You could therefore calculate the 
102 result once only and store it as a private data member of the class. 
103 It is optional whether you want to use this method, however, or put 
104 everything in <code>sigmaHat</code>.
105
106 <method name="double sigmaHat()">
107 is the key method for cross section calculations and returns a cross section
108 value, as further described below. It is called when also a preliminary set 
109 of incoming flavours has been picked, in addition to the kinematical ones 
110 already available for <code>sigmaKin</code>. Typically <code>sigmaHat</code> 
111 is called inside a loop over all allowed incoming flavour combinations,
112 stored in <code>id1</code> and <code>id2</code>, with fixed kinematics, 
113 as already illustrated above. The sum over the different flavour combinations 
114 provides the total cross section, while their relative size is used to make 
115 a selection of a specific incomimg state.
116 <br/>For a <ei>2 -> 1</ei> process, the returned value should be 
117 <ei>sigmaHat(sHat)</ei>, where <code>mH</code> (= <ei>mHat</ei>), 
118 <code>sH</code> (= <ei>sHat</ei>) and <code>sH2</code> (= <ei>sHat^2</ei>) 
119 are available to be used. 
120 <br/>For a <ei>2 -> 2</ei> process, instead
121 <ei>d(sigmaHat)/d(tHat)</ei> should be returned, based on
122 provided <code>mH, sH, sH2, tH, tH2, uH, uH2, m3, s3, m4, s4</code> and 
123 <code>pT2</code> values (<code>s3 = m3*m3</code> etc.). 
124 <br/>For a <ei>2 -> 3</ei> process, instead <ei>|M|^2</ei> should be 
125 returned, with normalization such that <ei>|M|^2 / (2 sHat)</ei> integrated 
126 over the three-body phase space gives the cross section. Here no standard 
127 set of variables exist. Instead the obvious ones, 
128 <code>mH, sH, m3, s3, m4, s4, m5, s5</code>, are complemented by the 
129 four-vectors <code>p3cm, p4cm, p5cm</code>, from which further invariants 
130 may be calculated. The four-vectors are defined in the cm frame of the 
131 subcollision, with incoming partons along the <ei>+-z</ei> axis.   
132 <br/>In either case, <ei>alpha_s</ei> and <ei>alpha_em</ei> have already 
133 been calculated, and are stored in <code>alpS</code> and <code>alpEM</code>. 
134 Also other standard variables may be used, like 
135 <code>CoupEW::sin2thetaW()</code>, and related flavour-dependent
136 vector and axial couplings in <code>CoupEW</code> and CKM combinations
137 in <code>VCKM</code>. 
138 <br/>In case some of the final-state particles are resonances, their 
139 squared masses have already been selected according to a Breit-Wigner
140 with a linearly running width <ei>Gamma(m) = Gamma(m_0) * m / m_0</ei>.
141 More precisely, the mass spectrum is weighted according to
142 <ei>w_BW(m^2) d(m^2)</ei>, where
143 <eq>
144 w_BW(m^2) = (1/pi) * (m * Gamma(m)) / ( (m^2 - m_0^2)^2 + (m * Gamma(m))^2 ) . 
145 </eq> 
146 If you would like to have another expression, the above weights are stored 
147 in <code>runBW3</code>, <code>runBW4</code> and <code>runBW5</code>, 
148 respectively. If you divide out one of these factors, you just remain with 
149 a phase space selection <ei>d(m^2)</ei> for this particle,
150 and can multiply on your desired shape factor instead. Unfortunately, the 
151 Monte Carlo efficiency will drop if your new mass distribution differs 
152 dramatically from the input one. Therefore it does make sense to adjust the 
153 database value of the width to be slightly (but not too much) broader 
154 than the distribution you have in mind. Also note that, already by default, 
155 the wings of the Breit-Wigner are oversampled (with a compensating lower 
156 internal weight) by partly sampling like <ei>(a + b/m^2 + c/m^4) d(m^2)</ei>,
157 where the last term is only used for <ei>gamma^*/Z^0</ei>.
158
159 <method name="void setIdColAcol()">
160 is called only once an initial state and a kinematical configuration has 
161 been picked. This routine must set the complete flavour information and 
162 the colour flow of the process. This may involve further random choices,
163 between different possible final-state flavours or between possible 
164 competing colour flows. Private data members of the class may be used to 
165 retain some information from the previous steps above.
166 <br/>When this routine is called the two incoming flavours have already 
167 been selected and are available in <code>id1</code> and <code>id2</code>, 
168 whereas the one, two or three outgoing ones either are fixed for a given 
169 process or can be determined from the instate (e.g. whether a <ei>W^+</ei> 
170 or <ei>W^-</ei> was produced).  There is also a standard method in 
171 <code>VCKM</code> to pick a final flavour from an initial one with CKM 
172 mixing. Once you have figured out the value of
173 <code>id3</code> and, the case being, <code>id4</code> and 
174 <code>id5</code>, you store these values permanently by a call
175 <code>setId( id1, id2, id3, id4, id5)</code>, where the last two may be 
176 omitted if irrelevant. 
177 <br/>Correspondingly, the colours are stored with
178 <code>setColAcol( col1, acol1, col2, acol2, col3, acol3, col4, acol4, 
179 col5, acol5)</code>, where the final ones may be omitted if irrelevant.
180 Les Houches style colour tags are used, but starting with number 1
181 (and later shifted by the currently requested offset). The 
182 input is grouped particle by particle, with the colour index before the 
183 anticolour one. You may need to select colour flow dynamically, depending 
184 on the kinematics, when several distinct possibilities exist. Trivial 
185 operations, like swapping colours and anticolours, can be done with 
186 existing methods.
187 <br/>When the <code>id3Mass()</code> and <code>id4Mass()</code>
188 methods have been used, the order of the outgoing particles may be 
189 inconsistent with the way the <ei>tHat</ei> and <ei>uHat</ei>
190 variables have been defined. A typical example would be a process like
191 <ei>q g -> q' W</ei> with <ei>tHat</ei> defined between incoming and 
192 outgoing quark, but where <code>id3Mass() = 24</code> and so the
193 process is to be stored as <ei>q g -> W q'</ei>. One should then put
194 the variable <code>swapTU = true</code> in <code>setIdColAcol()</code>
195 for each event where the <ei>tHat</ei> and <ei>uHat</ei> variables 
196 should be swapped before the event kinematics is reconstructed. This 
197 variable is automatically restored to <code>false</code> for each new 
198 event.
199
200 <method name="double weightDecayFlav( Event& process)">
201 is called to allow a reweighting of the simultaneous flavour choices of 
202 resonance decay products. Is currently only used for the 
203 <ei>q qbar -> gamma*/Z^0 gamma*/Z^0</ei> process, and will likely not
204 be of interest for you. 
205
206 <method name="double weightDecay( Event& process, int iResBeg, int iResEnd)">
207 is called when the basic process has one or several resonances, after each 
208 set of related resonances in <code>process[i]</code>,
209 <code>iResBeg</code> &lt;= <code>i </code> &lt;= <code>iResEnd</code>, 
210 has been allowed to decay. The calculated weight, to be normalized 
211 to the range between 0 and 1, is used to decide whether to accept the 
212 decay(s) or try for a new decay configuration. The base-class version of
213 this method returns unity, i.e. gives isotropic decays by default. 
214 This method may be called repeatedly for a single event. For instance, in 
215 <ei>q qbar -> H^0 Z^0</ei> with <ei>H^0 -> W^+ W^-</ei>, a first call 
216 would be made after the <ei>H^0</ei> and <ei>Z^0</ei> decays, and then 
217 depend only on the <ei>Z^0</ei> decay angles since the <ei>H^0</ei> 
218 decays isotropically. The second call would be after the <ei>W^+ W^-</ei> 
219 decays and then involve correlations between the four daughter fermions.
220
221 <method name="string name()">
222 returns the name of the process, as you want it to be shown in listings.
223
224 <method name="int code()"> 
225 returns an integer identifier of the process. This has no internal function, 
226 but is only intended as a service for the user to rapidly (and hopefully
227 uniquely) identify which process occured in a given event. Numbers below 
228 10000 are reserved for internal PYTHIA use. 
229
230 <method name="string inFlux()">
231 this string specifies the combinations of incoming partons that are 
232 allowed for the process under consideration, and thereby which incoming
233 flavours <code>id1</code> and <code>id2</code> the <code>sigmaHat()</code> 
234 calls will be looped over. It is always possible to pick a wider flavour 
235 selection than strictly required and then put to zero cross sections in 
236 the superfluous channels, but of course this may cost some extra execution 
237 time. Currently allowed options are:
238 <br/>* <code>gg</code>: two gluons. 
239 <br/>* <code>qg</code>: one (anti)quark and one gluon.
240 <br/>* <code>qq</code>: any combination of two quarks, two antiquarks or
241 a quark and an antiquark.
242 <br/>* <code>qqbarSame</code>: a quark and its antiquark;
243 this is a subset of the above <code>qq</code> option.
244 <br/>* <code>ff</code>: any combination of two fermions, two antifermions 
245 or a fermion and an antifermion; is the same as <code>qq</code> for 
246 hadron beams but also allows processes to work with lepton beams.
247 <br/>* <code>ffbarSame</code>: a fermion and its antifermion; is the 
248 same as <code>qqbarSame</code> for hadron beams but also allows processes 
249 to work with lepton beams.
250 <br/>* <code>ffbarChg</code>: a fermion and an antifermion that combine 
251 to give charge +-1.
252 <br/>* <code>fgm</code>: a fermion and a photon (gamma).
253 <br/>* <code>ggm</code>: a gluon and a photon.
254 <br/>* <code>gmgm</code>: two photons.
255
256 <method name="bool convert2mb()">
257 it is assumed that cross sections normally come in dimensions such that
258 they, when integrated over the relevant phase space, obtain the dimension
259 GeV^-2, and therefore need to be converted to mb. If the cross section 
260 is already encoded as mb then <code>convert2mb()</code> should be 
261 overloaded to instead return <code>false</code>.
262
263 <method name="int id3Mass(), int id4Mass(), int id5Mass()"> 
264 are the one, two or three final-state flavours, where masses are to be 
265 selected before the matrix elements are evaluated. Only the absolute value 
266 should be given. For massless particles, like gluons and photons, one need 
267 not give anything, i.e. one defaults to 0. The same goes for normal light 
268 quarks, where masses presumably are not implemented in the matrix elements.  
269 Later on, these quarks can still (automatically) obtain constituent masses, 
270 once a <ei>u</ei>, <ei>d</ei> or <ei>s</ei> flavour has been selected. 
271
272 <method name="int resonanceA(), int resonanceB()"> 
273 are the codes of up to two <ei>s</ei>-channel resonances contributing to 
274 the matrix elements. These are used by the program to improve the phase-space 
275 selection efficiency, by partly sampling according to the relevant 
276 Breit-Wigners. Massless resonances (the gluon and photon) need not be 
277 specified.
278
279 <method name="bool isSChannel()"> 
280 normally the choice of renormalization and factorization scales in 
281 <ei>2 -> 2</ei> and <ei>2 -> 3</ei> processes is based on the assumption 
282 that <ei>t</ei>- and <ei>u</ei>-channel exchanges dominates the 
283 cross section. In cases such as <ei>f fbar -> gamma* -> f' fbar'</ei> a 
284 <ei>2 -> 2</ei> process actually ought to be given scales as a 
285 <ei>2 -> 1</ei> one, in the sense that it proceeds entirely through 
286 an <ei>s</ei>-channel resonance. This can be achieved if you override the
287 default <code>false</code> to return <code>true</code>. See further the
288 page on <aloc href="CouplingsAndScales">couplings and scales</aloc>.
289
290 <method name="int idTchan1(), int idTchan2()">
291 the <ei>2 -> 3</ei> phase space selection machinery is rather primitive,
292 as already mentioned. The efficiency can be improved in processes that
293 proceed though <ei>t</ei>-channel exchanges, such as 
294 <ei>q qbar' -> H^0 q qbar'</ei> via <ei>Z^0 Z^0</ei> fusion, if the identity 
295 of the  <ei>t</ei>-channel-exchanged particles on the two side of the 
296 event are provided. Only the absolute value is of interest.
297
298 <method name="double tChanFracPow1(), double tChanFracPow2()">
299 in the above kind of <ei>2 -> 3</ei> phase-space selection, the
300 sampling of <ei>pT^2</ei> is done with one part flat, one part weighted
301 like <ei>1 / (pT^2 + m_R^2)</ei> and one part  like 
302 <ei>1 / (pT^2 + m_R^2)^2</ei>. The above values provide the relative
303 amount put in the latter two channels, respectively, with the first 
304 obtaining the rest. Thus the sum of <code>tChanFracPow1()</code> and
305 <code>tChanFracPow2()</code> must be below unity. The final results
306 should be independent of these numbers, but the Monte Carlo efficiency 
307 may be quite low for a bad choice. Here <ei>m_R</ei> is the mass of the 
308 exchanged resonance specified by <code>idTchan1()</code> or 
309 <code>idTchan2()</code>. Note that the order of the final-state 
310 listing is important in the above <ei>q qbar' -> H^0 q qbar'</ei> example, 
311 i.e. the <ei>H^0</ei> must be returned by <code>id3Mass()</code>,
312 since it is actually the <ei>pT^2</ei> of the latter two that are 
313 selected independently, with the first <ei>pT</ei> then fixed  
314 by transverse-momentum conservation.
315
316 <method name="useMirrorWeight()">
317 in <ei>2 -> 3</ei> processes the phase space selection used here
318 involves a twofold ambiguity basically corresponding to a flipping of 
319 the positions of last two outgoing particles. These are assumed equally 
320 likely by default, <code>false</code>, but for processes proceeding entirely 
321 through <ei>t</ei>-channel exchange the Monte Carlo efficiency can be 
322 improved by making a preselection based on the relative propagator
323 weights, <code>true</code>.  
324
325 <method name="int gmZmode()">
326 allows a possibility to override the global mode 
327 <aloc href="ElectroweakProcesses"><code>WeakZ0:gmZmode</code></aloc> 
328 for a specific process. The global mode normally is used to switch off 
329 parts of the <ei>gamma^*/Z^0</ei> propagator for test purposes. The
330 above local mode is useful for processes where a <ei>Z^0</ei> really is
331 that and nothing more, such as <ei>q qbar -> H^0 Z^0</ei>. The default
332 value -1 returned by <code>gmZmode()</code> ensures that the global
333 mode is used. 
334
335 <h3>Access to a process</h3> 
336
337 Once you have implemented a class, it is straightforward to make use of 
338 it in a run. Assume you have written a new class <code>MySigma</code>, 
339 which inherits from <code>Sigma1Process</code>, <code>Sigma2Process</code> 
340 or <code>Sigma3Process</code>, which in their turn inherit from 
341 <code>SigmaProcess</code>. You then create an instance of this class
342 and hand it in to a <code>pythia</code> object with 
343 <pre>
344       SigmaProcess* mySigma = new MySigma();
345       pythia.setSigmaPtr( mySigma); 
346 </pre>
347 If you have several processes you can repeat the procedure any number
348 of times. When <code>pythia.init(...)</code> is called these processes 
349 are initialized along with any internal processes you may have switched on, 
350 and treated in exactly the same manner. The  <code>pythia.next()</code> 
351 will therefore generate a mix of the different kinds of processes without 
352 distinction. See also the <aloc href="ProgramFlow">Program Flow</aloc> 
353 description.
354
355 <p/>
356 If the code should be of good quality and general usefulness, it would
357 be simple to include it as a permanently available process in the 
358 standard program distribution. The final step of that integration ought to 
359 be left for the PYTHIA authors, but here is a description of what is 
360 required.
361  
362 <p/>
363 A flag has to be defined, that allows the process to be switched on;
364 by default it should always be off. The name of the flag should be 
365 chosen of the type <code>model:process</code>. Here the 
366 <code>model</code> would be related to the general scenario considered,
367 e.g. <code>Compositeness</code>, while <code>process</code> would
368 specify instate and outstate, separated by a 2 (= to), e.g.
369 <code>ug2u*g</code>. 
370 When several processes are implemented and "belong together" it is 
371 also useful to define a <code>model:all</code> switch that affects
372 all the separate processes. 
373
374 <p/>
375 The flags should normally be stored in the <code>ProcessSelection.xml</code>
376 file or one of its daughters for a specific kind of processes. This is to 
377 make them easily found by users. You could create and use your own 
378 <code>.xml</code> file, so long as you then add that name to the 
379 list of files in the <code>Index.xml</code> file. (If not,
380 the flags would never be created and the program would not work.)  
381
382 <p/>
383 In the <code>ProcessContainer.c</code> file, the
384 <code>SetupContainers::init()</code> method needs to be expanded to
385 create instances of the processes switched on. This code is fairly
386 repetitive, and should be easy to copy and modify from the code 
387 already there. The basic structure is 
388 <br/>(i) check whether a process is requested by the user and, if so, 
389 <br/>(ii) create an instance of the matrix-element class, 
390 <br/>(iii)create a container for the matrix element and its associated 
391 phase-space handling, and 
392 <br>(iv) add the container to the existing process list.  
393
394 <p/>
395 Two minor variations are possible. One is that a set of related 
396 processes are lumped inside the the same initial check, i.e. are 
397 switched on all together. The second is that the matrix-element 
398 constructor may take arguments, as specified by you (see above). 
399 If so, the same basic matrix element may be recycled for a set of 
400 related processes, e.g. one for a composite <ei>u</ei> and one for 
401 a composite <ei>d</ei>. Obviously these variations may be combined.
402
403 </chapter>
404
405 <!-- Copyright (C) 2008 Torbjorn Sjostrand -->