- Update to pythia8140
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / xmldoc / EventAnalysis.xml
1 <chapter name="Event Analysis">
2
3 <h2>Event Analysis</h2>
4
5 <h3>Introduction</h3>
6
7 The routines in this section are intended to be used to analyze
8 event properties. As such they are not part of the main event
9 generation chain, but can be used in comparisons between Monte
10 Carlo events and real data. They are rather free-standing, but 
11 assume that input is provided in the PYTHIA 8 
12 <code>Event</code> format, and use a few basic facilities such 
13 as four-vectors.
14
15 <p/>
16 In addition to the methods presented here, there is also the 
17 possibility to make use of <aloc href="JetFinders">external 
18 jet finders </aloc>.
19
20 <h3>Sphericity</h3>
21
22 The standard sphericity tensor is
23 <eq>
24     S^{ab} = (sum_i p_i^a p_i^b) / (sum_i p_i^2)
25 </eq>
26 where the <ei>sum i</ei> runs over the particles in the event,
27 <ei>a, b = x, y, z,</ei> and <ei>p</ei> without such an index is 
28 the absolute size of the three-momentum . This tensor can be 
29 diagonalized to find eigenvalues and eigenvectors.
30
31 <p/>
32 The above tensor can be generalized by introducing a power 
33 <ei>r</ei>, such that
34 <eq>
35     S^{ab} = (sum_i p_i^a p_i^b p_i^{r-2}) / (sum_i p_i^r)
36 </eq>
37 In particular, <ei>r = 1</ei> gives a linear dependence on momenta 
38 and thus a collinear safe definition, unlike sphericity.
39
40 <p/>  
41 To do sphericity analyses you have to set up a <code>Sphericity</code>
42 instance, and then feed in events to it, one at a time. The results 
43 for the latest event are available as output from a few methods.
44
45 <method name="Sphericity::Sphericity(double power = 2., int select = 2)">
46 create a sphericity analysis object, where
47 <argument name="power" default="2."> 
48 is the power <ei>r</ei> defined above, i.e. 
49 <argoption value="2.">gives Spericity, and</argoption> 
50 <argoption value="1.">gives the linear form.</argoption>
51 </argument>
52 <argument name="select" default="2"> 
53 tells which particles are analyzed,
54 <argoption value="1">all final-state particles,</argoption>
55 <argoption value="2">all observable final-state particles, 
56 i.e. excluding neutrinos and other particles without strong or
57 electromagnetic interactions (the <code>isVisible()</code> 
58 particle method), and
59 </argoption>
60 <argoption value="3">only charged final-state particles.</argoption>
61 </argument>
62 </method>
63
64 <method name="bool Sphericity::analyze( const Event& event, 
65 ostream& os = cout)">
66 perform a sphericity analysis, where 
67 <argument name="event">is an object of the <code>Event</code> class, 
68 most likely the <code>pythia.event</code> one.
69 </argument>
70 <argument name="os" default="cout"> is the output stream for 
71 error messages. (The method does not rely on the <code>Info</code>
72 mchinery for error messages.)
73 </argument>
74 <br/>If the routine returns <code>false</code> the 
75 analysis failed, e.g. if too few particles are present to analyze.
76 </method>
77
78 <p/>
79 After the analysis has been performed, a few methods are available 
80 to return the result of the analysis of the latest event:
81
82 <method name="double Sphericity::sphericity()">
83 gives the sphericity (or equivalent if <ei>r</ei> is not 2),
84 </method>
85
86 <method name="double Sphericity::aplanarity()"> 
87 gives the aplanarity (with the same comment),
88 </method>
89
90 <method name="double Sphericity::eigenValue(int i)"> 
91 gives one of the three eigenvalues for <ei>i</ei> = 1, 2 or 3, in 
92 descending order,
93 </method>
94
95 <method name="Vec4 Sphericity::EventAxis(i)"> 
96 gives the matching normalized eigenvector, as a <code>Vec4</code> 
97 with vanishing time/energy component.
98 </method>
99
100 <method name="void Sphericity::list(ostream& os = cout)"> 
101 provides a listing of the above information.
102 </method>
103
104 <p/>
105 There is also one method that returns information accumulated for all
106 the events analyzed so far.
107
108 <method name="int Sphericity::nError()"> 
109 tells the number of times <code>analyze(...)</code> failed to analyze 
110 events, i.e. returned <code>false</code>.
111 </method>
112
113 <h3>Thrust</h3>
114
115 Thrust is obtained by varying the thrust axis so that the longitudinal
116 momentum component projected onto it is maximized, and thrust itself is 
117 then defined as the sum of absolute longitudinal momenta divided by
118 the sum of absolute momenta. The major axis is found correspondingly 
119 in the plane transverse to thrust, and the minor one is then defined 
120 to be transverse to both. Oblateness is the difference between the major 
121 and the minor values. 
122
123 <p/>
124 The calculation of thrust is more computer-time-intensive than e.g. 
125 linear sphericity, introduced above, and has no specific advantages except 
126 historical precedent. In the PYTHIA 6 implementation the search was 
127 speeded up at the price of then not being guaranteed to hit the absolute
128 maximum. The current implementation studies all possibilities, but at
129 the price of being slower, with time consumption for an event with
130 <ei>n</ei> particles growing like <ei>n^3</ei>.
131
132 <p/>  
133 To do thrust analyses you have to set up a <code>Thrust</code>
134 instance, and then feed in events to it, one at a time. The results 
135 for the latest event are available as output from a few methods.
136
137 <method name="Thrust::Thrust(int select = 2)">
138 create a thrust analysis object, where 
139 <argument name="select" default="2"> 
140 tells which particles are analyzed,
141 <argoption value="1">all final-state particles,</argoption>
142 <argoption value="2">all observable final-state particles, 
143 i.e. excluding neutrinos and other particles without strong or
144 electromagnetic interactions (the <code>isVisible()</code> 
145 particle method), and
146 </argoption>
147 <argoption value="3">only charged final-state particles.</argoption>
148 </argument>
149 </method>
150
151 <method name="bool Thrust::analyze( const Event& event, 
152 ostream& os = cout)">
153 perform a thrust analysis, where 
154 <argument name="event">is an object of the <code>Event</code> class, 
155 most likely the <code>pythia.event</code> one.
156 </argument>
157 <argument name="os" default="cout"> is the output stream for 
158 error messages. (The method does not rely on the <code>Info</code>
159 mchinery for error messages.)
160 </argument>
161 <br/>If the routine returns <code>false</code> the 
162 analysis failed, e.g. if too few particles are present to analyze.
163 </method>
164
165 <p/>
166 After the analysis has been performed, a few methods are available 
167 to return the result of the analysis of the latest event:
168
169 <method name="double Thrust::thrust()">
170 </method>
171 <methodmore name="double Thrust::tMajor()">
172 </methodmore>
173 <methodmore name="double Thrust::tMinor()">
174 </methodmore>
175 <methodmore name="double Thrust::oblateness()">
176 gives the thrust, major, minor and oblateness values, respectively, 
177 </methodmore>
178
179 <method name="Vec4 Thrust::eventAxis(int i)"> 
180 gives the matching normalized event-axis vectors, for <ei>i</ei> = 1, 2 or 3
181 corresponding to thrust, major or minor, as a <code>Vec4</code> with 
182 vanishing time/energy component.
183 </method>
184
185 <method name="void Thrust::list(ostream& os = cout)"> 
186 provides a listing of the above information.
187 </method>
188
189 <p/>
190 There is also one method that returns information accumulated for all
191 the events analyzed so far.
192
193 <method name="int Thrust::nError()"> 
194 tells the number of times <code>analyze(...)</code> failed to analyze 
195 events, i.e. returned <code>false</code>.
196 </method>
197
198 <h3>ClusterJet</h3>
199
200 <code>ClusterJet</code> (a.k.a. <code>LUCLUS</code> and 
201 <code>PYCLUS</code>) is a clustering algorithm of the type used for 
202 analyses of <ei>e^+e^-</ei> events, see the PYTHIA 6 manual. All 
203 visible particles in the events are clustered into jets. A few options 
204 are available for some well-known distance measures. Cutoff 
205 distances can either be given in terms of a scaled quadratic quantity 
206 like <ei>y = pT^2/E^2</ei> or an unscaled linear one like <ei>pT</ei>. 
207
208 <p/>  
209 To do jet finding analyses you have to set up a <code>ClusterJet</code>
210 instance, and then feed in events to it, one at a time. The results 
211 for the latest event are available as output from a few methods.
212
213 <method name="ClusterJet::ClusterJet(string measure = &quot;Lund&quot;, 
214 int select = 2, int massSet = 2, bool precluster = false, 
215 bool reassign = false)">
216 create a <code>ClusterJet</code> instance, where 
217 <argument name="measure" default="&quot;Lund&quot;">distance measure, 
218 to be provided as a character string (actually, only the first character 
219 is necessary)
220 <argoption value="&quot;Lund&quot;">the Lund <ei>pT</ei> distance,
221 </argoption>
222 <argoption value="&quot;JADE&quot;">the JADE mass distance, and
223 </argoption>
224 <argoption value="&quot;Durham&quot;">the Durham <ei>kT</ei> measure.
225 </argoption>
226 </argument>
227 <argument name="select" default="2"> 
228 tells which particles are analyzed,
229 <argoption value="1">all final-state particles,</argoption>
230 <argoption value="2">all observable final-state particles, 
231 i.e. excluding neutrinos and other particles without strong or
232 electromagnetic interactions (the <code>isVisible()</code> particle 
233 method), and
234 </argoption>
235 <argoption value="3">only charged final-state particles.</argoption>
236 </argument>
237 <argument name="massSet" default="2">masses assumed for the particles 
238 used in the analysis
239 <argoption value="0">all massless,</argoption>
240 <argoption value="1">photons are massless while all others are
241 assigned the <ei>pi+-</ei> mass, and
242 </argoption>
243 <argoption value="2">all given their correct masses.</argoption>
244 </argument>
245 <argument name="precluster" default="false">
246 perform or not a early preclustering step, where nearby particles
247 are lumped together so as to speed up the subsequent normal clustering.
248 </argument>
249 <argument name="reassign" default="false"> 
250 reassign all particles to the nearest jet each time after two jets
251 have been joined. 
252 </argument>
253 </method>
254
255 <method name="ClusterJet::analyze( const Event& event, double yScale, 
256 double pTscale, int nJetMin = 1, int nJetMax = 0, ostream& os = cout)"> 
257 performs a jet finding analysis, where 
258 <argument name="event">is an object of the <code>Event</code> class, 
259 most likely the <code>pythia.event</code> one.
260 </argument>
261 <argument name="yScale"> 
262 is the cutoff joining scale, below which jets are joined. Is given
263 in quadratic dimensionless quantities. Either <code>yScale</code>
264 or <code>pTscale</code> must be set nonvanishing, and the larger of 
265 the two dictates the actual value.
266 </argument>
267 <argument name="pTscale"> 
268 is the cutoff joining scale, below which jets are joined. Is given
269 in linear quantities, such as <ei>pT</ei> or <ei>m</ei> depending on
270 the measure used, but always in units of GeV. Either <code>yScale</code>
271 or <code>pTscale</code> must be set nonvanishing, and the larger of 
272 the two dictates the actual value.
273 </argument>
274 <argument name="nJetMin" default="1"> 
275 the minimum number of jets to be reconstructed. If used, it can override 
276 the <code>yScale</code> and <code>pTscale</code> values. 
277 </argument>
278 <argument name="nJetMax" default="0"> 
279 the maximum number of jets to be reconstructed. Is not used if below
280 <code>nJetMin</code>. If used, it can override the <code>yScale</code>
281 and <code>pTscale</code> values. Thus e.g. 
282 <code>nJetMin = nJetMax = 3</code> can be used to reconstruct exactly
283 3 jets.
284 </argument>
285 <argument name="os" default="cout"> is the output stream for 
286 error messages. (The method does not rely on the <code>Info</code>
287 mchinery for error messages.)
288 </argument>
289 <br/>If the routine returns <code>false</code> the analysis failed, 
290 e.g. because the number of particles was smaller than the minimum number
291 of jets requested.
292 </method>
293
294 <p/>
295 After the analysis has been performed, a few <code>ClusterJet</code> 
296 class methods are available to return the result of the analysis:
297
298 <method name="int ClusterJet::size()">
299 gives the number of jets found, with jets numbered 0 through 
300 <code>size() - 1</code>,
301 </method>
302
303 <method name="Vec4 ClusterJet::p(int i)">
304 gives a <code>Vec4</code> corresponding to the four-momentum defined by 
305 the sum of all the contributing particles to the <ei>i</ei>'th jet,
306 </method>
307
308 <method name="int ClusterJet::jetAssignment(int i)">
309 gives the index of the jet that the particle <ei>i</ei> of the event
310 record belongs to,
311 </method>
312
313 <method name="void ClusterJet::list(ostream& os = cout)">
314 provides a listing of the reconstructed jets.
315 </method>
316
317 <method name="int ClusterJet::distanceSize()">
318 the number of most recent clustering scales that have been stored
319 for readout with the next method. Normally this would be five, 
320 but less if fewer clustering steps occured.
321 </method>
322
323 <method name="double ClusterJet::distance(int i)">
324 clustering scales, with <code>distance(0)</code> being the most 
325 recent one, i.e. normally the highest, up to <code>distance(4)</code> 
326 being the fifth most recent. That is, with <ei>n</ei> being the final
327 number of jets, <code>ClusterJet::size()</code>, the scales at which
328 <ei>n+1</ei> jets become <ei>n</ei>, <ei>n+2</ei> become <ei>n+1</ei>,
329 and so on till <ei>n+5</ei> become <ei>n+4</ei>. Nonexisting clustering 
330 scales are returned as zero. The physical interpretation of a scale is 
331 as provided by the respective distance measure (Lund, JADE, Durham).
332 </method>
333
334 <p/>
335 There is also one method that returns information accumulated for all
336 the events analyzed so far.
337
338 <method name="int ClusterJet::nError()"> 
339 tells the number of times <code>analyze(...)</code> failed to analyze 
340 events, i.e. returned <code>false</code>.
341 </method>
342
343 <h3>CellJet</h3>
344
345 <code>CellJet</code> (a.k.a. <code>PYCELL</code>) is a simple cone jet 
346 finder in the UA1 spirit, see the PYTHIA 6 manual. It works in an 
347 <ei>(eta, phi, eT)</ei> space, where <ei>eta</ei> is pseudorapidity, 
348 <ei>phi</ei> azimuthal angle and <ei>eT</ei> transverse energy.
349 It will draw cones in <ei>R = sqrt(Delta-eta^2 + Delta-phi^2)</ei> 
350 around seed cells. If the total <ei>eT</ei> inside the cone exceeds 
351 the threshold, a jet is formed, and the cells are removed from further 
352 analysis. There are no split or merge procedures, so later-found jets 
353 may be missing some of the edge regions already used up by previous 
354 ones. Not all particles in the event are assigned to jets; leftovers 
355 may be viewed as belonging to beam remnants or the underlying event. 
356 It is not used by any experimental collaboration, but is closely 
357 related to the more recent and better theoretically motivated 
358 anti-<ei>kT</ei> algorithm <ref>Cac08</ref>.   
359
360 <p/>  
361 To do jet finding analyses you have to set up a <code>CellJet</code>
362 instance, and then feed in events to it, one at a time. Especially note 
363 that, if you want to use the options where energies are smeared in
364 order so emulate detector imperfections, you must hand in an external
365 random number generator, preferably the one residing in the 
366 <code>Pythia</code> class. The results for the latest event are 
367 available as output from a few methods.
368
369 <method name="CellJet::CellJet(double etaMax = 5., int nEta = 50, 
370 int nPhi = 32, int select = 2, int smear = 0, double resolution = 0.5, 
371 double upperCut = 2., double threshold = 0., Rndm* rndmPtr = 0)">
372 create a <code>CellJet</code> instance, where 
373 <argument name="etaMax" default="5."> 
374 the maximum +-pseudorapidity that the detector is assumed to cover.
375 </argument>
376 <argument name="nEta" default="50"> 
377 the number of equal-sized bins that the <ei>+-etaMax</ei> range 
378 is assumed to be divided into.
379 </argument>
380 <argument name="nPhi" default="32"> 
381 the number of equal-sized bins that the <ei>phi</ei> range 
382 <ei>+-pi</ei> is assumed to be divided into. 
383 </argument>
384 <argument name="select" default="2"> 
385 tells which particles are analyzed,
386 <argoption value="1">all final-state particles,</argoption>
387 <argoption value="2">all observable final-state particles, 
388 i.e. excluding neutrinos and other particles without strong or
389 electromagnetic interactions (the <code>isVisible()</code> particle 
390 method), 
391 and</argoption>
392 <argoption value="3">only charged final-state particles.</argoption>
393 </argument>
394 <argument name="smear" default="0">
395 strategy to smear the actual <ei>eT</ei> bin by bin, 
396 <argoption value="0">no smearing,</argoption>
397 <argoption value="1">smear the <ei>eT</ei> according to a Gaussian 
398 with width <ei>resolution * sqrt(eT)</ei>, with the Gaussian truncated 
399 at 0 and <ei>upperCut * eT</ei>,</argoption>
400 <argoption value="2">smear the <ei>e = eT * cosh(eta)</ei> according 
401 to a Gaussian with width <ei>resolution * sqrt(e)</ei>, with the 
402 Gaussian truncated at 0 and <ei>upperCut * e</ei>.</argoption>
403 </argument>
404 <argument name="resolution" default="0.5">
405 see above.
406 </argument>
407 <argument name="upperCut" default="2.">
408 see above.
409 </argument>
410 <argument name="threshold" default="0 GeV">
411 completely neglect all bins with an <ei>eT &lt; threshold</ei>.
412 </argument>
413 <argument name="rndmPtr" default="0">
414 the random-number generator used to select the smearing described
415 above. Must be handed in for smearing to be possible. If your 
416 <code>Pythia</code> class instance is named <code>pythia</code>,
417 then <code>&pythia.rndm</code> would be the logical choice.
418 </argument>
419 </class>
420
421 <method name="bool CellJet::analyze( const Event& event, 
422 double eTjetMin = 20., double coneRadius = 0.7, double eTseed = 1.5, 
423 ostream& os = cout)"> 
424 performs a jet finding analysis, where 
425 <argument name="event">is an object of the <code>Event</code> class, 
426 most likely the <code>pythia.event</code> one.
427 </argument>
428 <argument name="eTjetMin" default="20. GeV"> 
429 is the minimum transverse energy inside a cone for this to be 
430 accepted as a jet.
431 </argument>
432 <argument name="coneRadius" default="0.7"> 
433  is the size of the cone in <ei>(eta, phi)</ei> space drawn around 
434 the geometric center of the jet.
435 </argument>
436 <argument name="eTseed" default="1.5 GeV"> 
437 the mimimum <ei>eT</ei> in a cell for this to be acceptable as 
438 the trial center of a jet. 
439 </argument>
440 <argument name="os" default="cout"> is the output stream for 
441 error messages. (The method does not rely on the <code>Info</code>
442 mchinery for error messages.)
443 </argument>
444 <br/>If the routine returns <code>false</code> the analysis failed, 
445 but currently this is not foreseen ever to happen.
446 </method>
447
448 <p/>
449 After the analysis has been performed, a few <code>CellJet</code> 
450 class methods are available to return the result of the analysis:
451
452 <method name="int CellJet::size()">
453 gives the number of jets found, with jets numbered 0 through 
454 <code>size() - 1</code>,
455 </method>
456
457 <method name="double CellJet::eT(i)">
458 gives the <ei>eT</ei> of the <ei>i</ei>'th jet, where jets have been
459 ordered with decreasing <ei>eT</ei> values,
460 </method>
461
462 <method name="double CellJet::etaCenter(int i)">
463 </method>
464 <methodmore name="double CellJet::phiCenter(int i)">
465 gives the <ei>eta</ei> and <ei>phi</ei> coordinates of the geometrical 
466 center of the <ei>i</ei>'th jet,
467 </methodmore>
468
469 <method name="double CellJet::etaWeighted(int i)">
470 </method>
471 <methodmore name="double CellJet::phiWeighted(int i)">
472 gives the <ei>eta</ei> and <ei>phi</ei> coordinates of the 
473 <ei>eT</ei>-weighted center of the <ei>i</ei>'th jet,
474 </methodmore>
475
476 <method name="int CellJet::multiplicity(int i)">
477 gives the number of particles clustered into the <ei>i</ei>'th jet,
478 </method>
479
480 <method name="Vec4 CellJet::pMassless(int i)">
481 gives a <code>Vec4</code> corresponding to the four-momentum defined 
482 by the <ei>eT</ei> and the weighted center of the <ei>i</ei>'th jet,
483 </method>
484
485 <method name="Vec4 CellJet::pMassive(int i)">
486 gives a <code>Vec4</code> corresponding to the four-momentum defined by 
487 the sum of all the contributing cells to the <ei>i</ei>'th jet, where 
488 each cell contributes a four-momentum as if all the <ei>eT</ei> is 
489 deposited in the center of the cell,
490 </method>
491
492 <method name="double CellJet::m(int i)">
493 gives the invariant mass of the <ei>i</ei>'th jet, defined by the 
494 <code>pMassive</code> above,
495 </method>
496
497 <method name="void CellJet::list()">
498 provides a listing of the above information (except <code>pMassless</code>, 
499 for reasons of space).
500 </method>
501
502 <p/>
503 There is also one method that returns information accumulated for all
504 the events analyzed so far.
505 <method name="int CellJet::nError()"> 
506 tells the number of times <code>analyze(...)</code> failed to analyze 
507 events, i.e. returned <code>false</code>.
508 </method>
509
510 </chapter>
511
512 <!-- Copyright (C) 2010 Torbjorn Sjostrand -->