1 <chapter name="Semi-Internal Processes">
3 <h2>Semi-Internal Processes</h2>
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>.
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.
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.
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>.
43 <h3>The Cross Section Class</h3>
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.
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>
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.
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.
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,
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.
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>.
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
144 w_BW(m^2) = (1/pi) * (m * Gamma(m)) / ( (m^2 - m_0^2)^2 + (m * Gamma(m))^2 ) .
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>.
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
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
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.
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> <= <code>i </code> <= <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.
221 <method name="string name()">
222 returns the name of the process, as you want it to be shown in listings.
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.
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
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.
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>.
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.
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
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>.
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.
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.
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>.
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
335 <h3>Access to a process</h3>
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
344 SigmaProcess* mySigma = new MySigma();
345 pythia.setSigmaPtr( mySigma);
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>
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
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.
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.
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.)
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.
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.
405 <!-- Copyright (C) 2008 Torbjorn Sjostrand -->