1 <chapter name="Event Analysis">
3 <h2>Event Analysis</h2>
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
16 In addition to the methods presented here, there is also the
17 possibility to make use of <aloc href="JetFinders">external
22 The standard sphericity tensor is
24 S^{ab} = (sum_i p_i^a p_i^b) / (sum_i p_i^2)
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.
32 The above tensor can be generalized by introducing a power
35 S^{ab} = (sum_i p_i^a p_i^b p_i^{r-2}) / (sum_i p_i^r)
37 In particular, <ei>r = 1</ei> gives a linear dependence on momenta
38 and thus a collinear safe definition, unlike sphericity.
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.
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>
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>
60 <argoption value="3">only charged final-state particles.</argoption>
64 <method name="bool Sphericity::analyze( const Event& event,
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.
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.)
74 <br/>If the routine returns <code>false</code> the
75 analysis failed, e.g. if too few particles are present to analyze.
79 After the analysis has been performed, a few methods are available
80 to return the result of the analysis of the latest event:
82 <method name="double Sphericity::sphericity()">
83 gives the sphericity (or equivalent if <ei>r</ei> is not 2),
86 <method name="double Sphericity::aplanarity()">
87 gives the aplanarity (with the same comment),
90 <method name="double Sphericity::eigenValue(int i)">
91 gives one of the three eigenvalues for <ei>i</ei> = 1, 2 or 3, in
95 <method name="Vec4 Sphericity::eventAxis(i)">
96 gives the matching normalized eigenvector, as a <code>Vec4</code>
97 with vanishing time/energy component.
100 <method name="void Sphericity::list(ostream& os = cout)">
101 provides a listing of the above information.
105 There is also one method that returns information accumulated for all
106 the events analyzed so far.
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>.
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.
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>.
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.
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
147 <argoption value="3">only charged final-state particles.</argoption>
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.
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.)
161 <br/>If the routine returns <code>false</code> the
162 analysis failed, e.g. if too few particles are present to analyze.
166 After the analysis has been performed, a few methods are available
167 to return the result of the analysis of the latest event:
169 <method name="double Thrust::thrust()">
171 <methodmore name="double Thrust::tMajor()">
173 <methodmore name="double Thrust::tMinor()">
175 <methodmore name="double Thrust::oblateness()">
176 gives the thrust, major, minor and oblateness values, respectively,
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.
185 <method name="void Thrust::list(ostream& os = cout)">
186 provides a listing of the above information.
190 There is also one method that returns information accumulated for all
191 the events analyzed so far.
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>.
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>.
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.
213 <method name="ClusterJet::ClusterJet(string measure = "Lund",
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=""Lund"">distance measure,
218 to be provided as a character string (actually, only the first character
220 <argoption value=""Lund"">the Lund <ei>pT</ei> distance,
222 <argoption value=""JADE"">the JADE mass distance, and
224 <argoption value=""Durham"">the Durham <ei>kT</ei> measure.
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
235 <argoption value="3">only charged final-state particles.</argoption>
237 <argument name="massSet" default="2">masses assumed for the particles
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
243 <argoption value="2">all given their correct masses.</argoption>
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.
249 <argument name="reassign" default="false">
250 reassign all particles to the nearest jet each time after two jets
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.
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.
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.
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.
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
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.)
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
295 After the analysis has been performed, a few <code>ClusterJet</code>
296 class methods are available to return the result of the analysis:
298 <method name="int ClusterJet::size()">
299 gives the number of jets found, with jets numbered 0 through
300 <code>size() - 1</code>,
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,
308 <method name="int ClusterJet::jetAssignment(int i)">
309 gives the index of the jet that the particle <ei>i</ei> of the event
313 <method name="void ClusterJet::list(ostream& os = cout)">
314 provides a listing of the reconstructed jets.
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.
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).
335 There is also one method that returns information accumulated for all
336 the events analyzed so far.
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>.
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>.
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.
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.
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.
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.
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
392 <argoption value="3">only charged final-state particles.</argoption>
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>
404 <argument name="resolution" default="0.5">
407 <argument name="upperCut" default="2.">
410 <argument name="threshold" default="0 GeV">
411 completely neglect all bins with an <ei>eT < threshold</ei>.
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.
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.
428 <argument name="eTjetMin" default="20. GeV">
429 is the minimum transverse energy inside a cone for this to be
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.
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.
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.)
444 <br/>If the routine returns <code>false</code> the analysis failed,
445 but currently this is not foreseen ever to happen.
449 After the analysis has been performed, a few <code>CellJet</code>
450 class methods are available to return the result of the analysis:
452 <method name="int CellJet::size()">
453 gives the number of jets found, with jets numbered 0 through
454 <code>size() - 1</code>,
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,
462 <method name="double CellJet::etaCenter(int i)">
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,
469 <method name="double CellJet::etaWeighted(int i)">
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,
476 <method name="int CellJet::multiplicity(int i)">
477 gives the number of particles clustered into the <ei>i</ei>'th jet,
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,
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,
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,
497 <method name="void CellJet::list()">
498 provides a listing of the above information (except <code>pMassless</code>,
499 for reasons of space).
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>.
512 <!-- Copyright (C) 2010 Torbjorn Sjostrand -->