source: HiSusy/trunk/Pythia8/pythia8170/htmldoc/EventAnalysis.html @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES
File size: 40.6 KB
Line 
1<html>
2<head>
3<title>Event Analysis</title>
4<link rel="stylesheet" type="text/css" href="pythia.css"/>
5<link rel="shortcut icon" href="pythia32.gif"/>
6</head>
7<body>
8
9<h2>Event Analysis</h2>
10
11<h3>Introduction</h3>
12
13The routines in this section are intended to be used to analyze
14event properties. As such they are not part of the main event
15generation chain, but can be used in comparisons between Monte
16Carlo events and real data. They are rather free-standing, but
17assume that input is provided in the PYTHIA 8
18<code>Event</code> format, and use a few basic facilities such
19as four-vectors. Their ordering is mainly by history; for current
20LHC applications the final one, <code>SlowJet</code>, is of
21special interest.
22
23<p/>
24In addition to the methods presented here, there is also the
25possibility to make use of <a href="JetFinders.html" target="page">external
26jet finders </a>.
27
28<h3>Sphericity</h3>
29
30The standard sphericity tensor is
31<br/><i>
32    S^{ab} = (sum_i p_i^a p_i^b) / (sum_i p_i^2)
33</i><br/>
34where the <i>sum i</i> runs over the particles in the event,
35<i>a, b = x, y, z,</i> and <i>p</i> without such an index is
36the absolute size of the three-momentum . This tensor can be
37diagonalized to find eigenvalues and eigenvectors.
38
39<p/>
40The above tensor can be generalized by introducing a power
41<i>r</i>, such that
42<br/><i>
43    S^{ab} = (sum_i p_i^a p_i^b p_i^{r-2}) / (sum_i p_i^r)
44</i><br/>
45In particular, <i>r = 1</i> gives a linear dependence on momenta
46and thus a collinear safe definition, unlike sphericity.
47
48<p/> 
49To do sphericity analyses you have to set up a <code>Sphericity</code>
50instance, and then feed in events to it, one at a time. The results
51for the latest event are available as output from a few methods.
52
53<a name="method1"></a>
54<p/><strong>Sphericity::Sphericity(double power = 2., int select = 2) &nbsp;</strong> <br/>
55create a sphericity analysis object, where
56<br/><code>argument</code><strong> power </strong> (<code>default = <strong>2.</strong></code>) : 
57is the power <i>r</i> defined above, i.e.
58<br/><code>argumentoption </code><strong> 2.</strong> : gives Spericity, and   
59<br/><code>argumentoption </code><strong> 1.</strong> : gives the linear form. 
60 
61<br/><code>argument</code><strong> select </strong> (<code>default = <strong>2</strong></code>) : 
62tells which particles are analyzed,
63<br/><code>argumentoption </code><strong> 1</strong> : all final-state particles, 
64<br/><code>argumentoption </code><strong> 2</strong> : all observable final-state particles,
65i.e. excluding neutrinos and other particles without strong or
66electromagnetic interactions (the <code>isVisible()</code> 
67particle method), and
68 
69<br/><code>argumentoption </code><strong> 3</strong> : only charged final-state particles. 
70 
71 
72
73<a name="method2"></a>
74<p/><strong>bool Sphericity::analyze( const Event& event, ostream& os = cout) &nbsp;</strong> <br/>
75perform a sphericity analysis, where
76<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
77most likely the <code>pythia.event</code> one.
78 
79<br/><code>argument</code><strong> os </strong> (<code>default = <strong>cout</strong></code>) :  is the output stream for
80error messages. (The method does not rely on the <code>Info</code>
81mchinery for error messages.)
82 
83<br/>If the routine returns <code>false</code> the
84analysis failed, e.g. if too few particles are present to analyze.
85 
86
87<p/>
88After the analysis has been performed, a few methods are available
89to return the result of the analysis of the latest event:
90
91<a name="method3"></a>
92<p/><strong>double Sphericity::sphericity() &nbsp;</strong> <br/>
93gives the sphericity (or equivalent if <i>r</i> is not 2),
94 
95
96<a name="method4"></a>
97<p/><strong>double Sphericity::aplanarity() &nbsp;</strong> <br/>
98gives the aplanarity (with the same comment),
99 
100
101<a name="method5"></a>
102<p/><strong>double Sphericity::eigenValue(int i) &nbsp;</strong> <br/>
103gives one of the three eigenvalues for <i>i</i> = 1, 2 or 3, in
104descending order,
105 
106
107<a name="method6"></a>
108<p/><strong>Vec4 Sphericity::eventAxis(i) &nbsp;</strong> <br/>
109gives the matching normalized eigenvector, as a <code>Vec4</code> 
110with vanishing time/energy component.
111 
112
113<a name="method7"></a>
114<p/><strong>void Sphericity::list(ostream& os = cout) &nbsp;</strong> <br/>
115provides a listing of the above information.
116 
117
118<p/>
119There is also one method that returns information accumulated for all
120the events analyzed so far.
121
122<a name="method8"></a>
123<p/><strong>int Sphericity::nError() &nbsp;</strong> <br/>
124tells the number of times <code>analyze(...)</code> failed to analyze
125events, i.e. returned <code>false</code>.
126 
127
128<h3>Thrust</h3>
129
130Thrust is obtained by varying the thrust axis so that the longitudinal
131momentum component projected onto it is maximized, and thrust itself is
132then defined as the sum of absolute longitudinal momenta divided by
133the sum of absolute momenta. The major axis is found correspondingly
134in the plane transverse to thrust, and the minor one is then defined
135to be transverse to both. Oblateness is the difference between the major
136and the minor values.
137
138<p/>
139The calculation of thrust is more computer-time-intensive than e.g.
140linear sphericity, introduced above, and has no specific advantages except
141historical precedent. In the PYTHIA 6 implementation the search was
142speeded up at the price of then not being guaranteed to hit the absolute
143maximum. The current implementation studies all possibilities, but at
144the price of being slower, with time consumption for an event with
145<i>n</i> particles growing like <i>n^3</i>.
146
147<p/> 
148To do thrust analyses you have to set up a <code>Thrust</code>
149instance, and then feed in events to it, one at a time. The results
150for the latest event are available as output from a few methods.
151
152<a name="method9"></a>
153<p/><strong>Thrust::Thrust(int select = 2) &nbsp;</strong> <br/>
154create a thrust analysis object, where
155<br/><code>argument</code><strong> select </strong> (<code>default = <strong>2</strong></code>) : 
156tells which particles are analyzed,
157<br/><code>argumentoption </code><strong> 1</strong> : all final-state particles, 
158<br/><code>argumentoption </code><strong> 2</strong> : all observable final-state particles,
159i.e. excluding neutrinos and other particles without strong or
160electromagnetic interactions (the <code>isVisible()</code> 
161particle method), and
162 
163<br/><code>argumentoption </code><strong> 3</strong> : only charged final-state particles. 
164 
165 
166
167<a name="method10"></a>
168<p/><strong>bool Thrust::analyze( const Event& event, ostream& os = cout) &nbsp;</strong> <br/>
169perform a thrust analysis, where
170<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
171most likely the <code>pythia.event</code> one.
172 
173<br/><code>argument</code><strong> os </strong> (<code>default = <strong>cout</strong></code>) :  is the output stream for
174error messages. (The method does not rely on the <code>Info</code>
175mchinery for error messages.)
176 
177<br/>If the routine returns <code>false</code> the
178analysis failed, e.g. if too few particles are present to analyze.
179 
180
181<p/>
182After the analysis has been performed, a few methods are available
183to return the result of the analysis of the latest event:
184
185<a name="method11"></a>
186<p/><strong>double Thrust::thrust() &nbsp;</strong> <br/>
187 
188<strong>double Thrust::tMajor() &nbsp;</strong> <br/>
189 
190<strong>double Thrust::tMinor() &nbsp;</strong> <br/>
191 
192<strong>double Thrust::oblateness() &nbsp;</strong> <br/>
193gives the thrust, major, minor and oblateness values, respectively,
194 
195
196<a name="method12"></a>
197<p/><strong>Vec4 Thrust::eventAxis(int i) &nbsp;</strong> <br/>
198gives the matching normalized event-axis vectors, for <i>i</i> = 1, 2 or 3
199corresponding to thrust, major or minor, as a <code>Vec4</code> with
200vanishing time/energy component.
201 
202
203<a name="method13"></a>
204<p/><strong>void Thrust::list(ostream& os = cout) &nbsp;</strong> <br/>
205provides a listing of the above information.
206 
207
208<p/>
209There is also one method that returns information accumulated for all
210the events analyzed so far.
211
212<a name="method14"></a>
213<p/><strong>int Thrust::nError() &nbsp;</strong> <br/>
214tells the number of times <code>analyze(...)</code> failed to analyze
215events, i.e. returned <code>false</code>.
216 
217
218<h3>ClusterJet</h3>
219
220<code>ClusterJet</code> (a.k.a. <code>LUCLUS</code> and
221<code>PYCLUS</code>) is a clustering algorithm of the type used for
222analyses of <i>e^+e^-</i> events, see the PYTHIA 6 manual. All
223visible particles in the events are clustered into jets. A few options
224are available for some well-known distance measures. Cutoff
225distances can either be given in terms of a scaled quadratic quantity
226like <i>y = pT^2/E^2</i> or an unscaled linear one like <i>pT</i>.
227
228<p/> 
229Note that we have deliberately chosen not to include the <i>e^+e^-</i> 
230equivalents of the Cambridge/Aachen and anti-<i>kRT</i> algorithms.
231These tend to be good at clustering the densely populated (in angle)
232cores of jets, but less successful for the sparsely populated transverse
233regions, where many jets may come to consist of a single low-momentum
234particle. In hadron collisions such jets could easily be disregarded,
235while in <i>e^+e^-</i> annihilation all particles derive back from the
236hard process.
237
238<p/> 
239To do jet finding analyses you have to set up a <code>ClusterJet</code>
240instance, and then feed in events to it, one at a time. The results
241for the latest event are available as output from a few methods.
242
243<a name="method15"></a>
244<p/><strong>ClusterJet::ClusterJet(string measure = &quot;Lund&quot;, int select = 2, int massSet = 2, bool precluster = false, bool reassign = false) &nbsp;</strong> <br/>
245create a <code>ClusterJet</code> instance, where
246<br/><code>argument</code><strong> measure </strong> (<code>default = <strong>&quot;Lund&quot;</strong></code>) : distance measure,
247to be provided as a character string (actually, only the first character
248is necessary)
249<br/><code>argumentoption </code><strong> &quot;Lund&quot;</strong> : the Lund <i>pT</i> distance,
250 
251<br/><code>argumentoption </code><strong> &quot;JADE&quot;</strong> : the JADE mass distance, and
252 
253<br/><code>argumentoption </code><strong> &quot;Durham&quot;</strong> : the Durham <i>kT</i> measure.
254 
255 
256<br/><code>argument</code><strong> select </strong> (<code>default = <strong>2</strong></code>) : 
257tells which particles are analyzed,
258<br/><code>argumentoption </code><strong> 1</strong> : all final-state particles, 
259<br/><code>argumentoption </code><strong> 2</strong> : all observable final-state particles,
260i.e. excluding neutrinos and other particles without strong or
261electromagnetic interactions (the <code>isVisible()</code> particle
262method), and
263 
264<br/><code>argumentoption </code><strong> 3</strong> : only charged final-state particles. 
265 
266<br/><code>argument</code><strong> massSet </strong> (<code>default = <strong>2</strong></code>) : masses assumed for the particles
267used in the analysis
268<br/><code>argumentoption </code><strong> 0</strong> : all massless, 
269<br/><code>argumentoption </code><strong> 1</strong> : photons are massless while all others are
270assigned the <i>pi+-</i> mass, and
271 
272<br/><code>argumentoption </code><strong> 2</strong> : all given their correct masses. 
273 
274<br/><code>argument</code><strong> precluster </strong> (<code>default = <strong>false</strong></code>) :
275perform or not a early preclustering step, where nearby particles
276are lumped together so as to speed up the subsequent normal clustering.
277 
278<br/><code>argument</code><strong> reassign </strong> (<code>default = <strong>false</strong></code>) : 
279reassign all particles to the nearest jet each time after two jets
280have been joined.
281 
282 
283
284<a name="method16"></a>
285<p/><strong>ClusterJet::analyze( const Event& event, double yScale, double pTscale, int nJetMin = 1, int nJetMax = 0, ostream& os = cout) &nbsp;</strong> <br/>
286performs a jet finding analysis, where
287<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
288most likely the <code>pythia.event</code> one.
289 
290<br/><code>argument</code><strong> yScale </strong>  : 
291is the cutoff joining scale, below which jets are joined. Is given
292in quadratic dimensionless quantities. Either <code>yScale</code>
293or <code>pTscale</code> must be set nonvanishing, and the larger of
294the two dictates the actual value.
295 
296<br/><code>argument</code><strong> pTscale </strong>  : 
297is the cutoff joining scale, below which jets are joined. Is given
298in linear quantities, such as <i>pT</i> or <i>m</i> depending on
299the measure used, but always in units of GeV. Either <code>yScale</code>
300or <code>pTscale</code> must be set nonvanishing, and the larger of
301the two dictates the actual value.
302 
303<br/><code>argument</code><strong> nJetMin </strong> (<code>default = <strong>1</strong></code>) : 
304the minimum number of jets to be reconstructed. If used, it can override
305the <code>yScale</code> and <code>pTscale</code> values.
306 
307<br/><code>argument</code><strong> nJetMax </strong> (<code>default = <strong>0</strong></code>) : 
308the maximum number of jets to be reconstructed. Is not used if below
309<code>nJetMin</code>. If used, it can override the <code>yScale</code>
310and <code>pTscale</code> values. Thus e.g.
311<code>nJetMin = nJetMax = 3</code> can be used to reconstruct exactly
3123 jets.
313 
314<br/><code>argument</code><strong> os </strong> (<code>default = <strong>cout</strong></code>) :  is the output stream for
315error messages. (The method does not rely on the <code>Info</code>
316mchinery for error messages.)
317 
318<br/>If the routine returns <code>false</code> the analysis failed,
319e.g. because the number of particles was smaller than the minimum number
320of jets requested.
321 
322
323<p/>
324After the analysis has been performed, a few <code>ClusterJet</code> 
325class methods are available to return the result of the analysis:
326
327<a name="method17"></a>
328<p/><strong>int ClusterJet::size() &nbsp;</strong> <br/>
329gives the number of jets found, with jets numbered 0 through
330<code>size() - 1</code>.
331 
332
333<a name="method18"></a>
334<p/><strong>Vec4 ClusterJet::p(int i) &nbsp;</strong> <br/>
335gives a <code>Vec4</code> corresponding to the four-momentum defined by
336the sum of all the contributing particles to the <i>i</i>'th jet.
337 
338
339<a name="method19"></a>
340<p/><strong>int ClusterJet::mult(int i) &nbsp;</strong> <br/>
341the number of particles that have been clustered into the <i>i</i>'th jet.
342 
343
344<a name="method20"></a>
345<p/><strong>int ClusterJet::jetAssignment(int i) &nbsp;</strong> <br/>
346gives the index of the jet that the particle <i>i</i> of the event
347record belongs to,
348 
349
350<a name="method21"></a>
351<p/><strong>void ClusterJet::list(ostream& os = cout) &nbsp;</strong> <br/>
352provides a listing of the reconstructed jets.
353 
354
355<a name="method22"></a>
356<p/><strong>int ClusterJet::distanceSize() &nbsp;</strong> <br/>
357the number of most recent clustering scales that have been stored
358for readout with the next method. Normally this would be five,
359but less if fewer clustering steps occured.
360 
361
362<a name="method23"></a>
363<p/><strong>double ClusterJet::distance(int i) &nbsp;</strong> <br/>
364clustering scales, with <code>distance(0)</code> being the most
365recent one, i.e. normally the highest, up to <code>distance(4)</code> 
366being the fifth most recent. That is, with <i>n</i> being the final
367number of jets, <code>ClusterJet::size()</code>, the scales at which
368<i>n+1</i> jets become <i>n</i>, <i>n+2</i> become <i>n+1</i>,
369and so on till <i>n+5</i> become <i>n+4</i>. Nonexisting clustering
370scales are returned as zero. The physical interpretation of a scale is
371as provided by the respective distance measure (Lund, JADE, Durham).
372 
373
374<p/>
375There is also one method that returns information accumulated for all
376the events analyzed so far.
377
378<a name="method24"></a>
379<p/><strong>int ClusterJet::nError() &nbsp;</strong> <br/>
380tells the number of times <code>analyze(...)</code> failed to analyze
381events, i.e. returned <code>false</code>.
382 
383
384<h3>CellJet</h3>
385
386<code>CellJet</code> (a.k.a. <code>PYCELL</code>) is a simple cone jet
387finder in the UA1 spirit, see the PYTHIA 6 manual. It works in an
388<i>(eta, phi, eT)</i> space, where <i>eta</i> is pseudorapidity,
389<i>phi</i> azimuthal angle and <i>eT</i> transverse energy.
390It will draw cones in <i>R = sqrt(Delta-eta^2 + Delta-phi^2)</i> 
391around seed cells. If the total <i>eT</i> inside the cone exceeds
392the threshold, a jet is formed, and the cells are removed from further
393analysis. There are no split or merge procedures, so later-found jets
394may be missing some of the edge regions already used up by previous
395ones. Not all particles in the event are assigned to jets; leftovers
396may be viewed as belonging to beam remnants or the underlying event.
397It is not used by any experimental collaboration, but is closely
398related to the more recent and better theoretically motivated
399anti-<i>kT</i> algorithm [<a href="Bibliography.html" target="page">Cac08</a>].   
400
401<p/> 
402To do jet finding analyses you have to set up a <code>CellJet</code>
403instance, and then feed in events to it, one at a time. Especially note
404that, if you want to use the options where energies are smeared in
405order so emulate detector imperfections, you must hand in an external
406random number generator, preferably the one residing in the
407<code>Pythia</code> class. The results for the latest event are
408available as output from a few methods.
409
410<a name="method25"></a>
411<p/><strong>CellJet::CellJet(double etaMax = 5., int nEta = 50, int nPhi = 32, int select = 2, int smear = 0, double resolution = 0.5, double upperCut = 2., double threshold = 0., Rndm* rndmPtr = 0) &nbsp;</strong> <br/>
412create a <code>CellJet</code> instance, where
413<br/><code>argument</code><strong> etaMax </strong> (<code>default = <strong>5.</strong></code>) : 
414the maximum +-pseudorapidity that the detector is assumed to cover.
415 
416<br/><code>argument</code><strong> nEta </strong> (<code>default = <strong>50</strong></code>) : 
417the number of equal-sized bins that the <i>+-etaMax</i> range
418is assumed to be divided into.
419 
420<br/><code>argument</code><strong> nPhi </strong> (<code>default = <strong>32</strong></code>) : 
421the number of equal-sized bins that the <i>phi</i> range
422<i>+-pi</i> is assumed to be divided into.
423 
424<br/><code>argument</code><strong> select </strong> (<code>default = <strong>2</strong></code>) : 
425tells which particles are analyzed,
426<br/><code>argumentoption </code><strong> 1</strong> : all final-state particles, 
427<br/><code>argumentoption </code><strong> 2</strong> : all observable final-state particles,
428i.e. excluding neutrinos and other particles without strong or
429electromagnetic interactions (the <code>isVisible()</code> particle
430method),
431and 
432<br/><code>argumentoption </code><strong> 3</strong> : only charged final-state particles. 
433 
434<br/><code>argument</code><strong> smear </strong> (<code>default = <strong>0</strong></code>) :
435strategy to smear the actual <i>eT</i> bin by bin,
436<br/><code>argumentoption </code><strong> 0</strong> : no smearing, 
437<br/><code>argumentoption </code><strong> 1</strong> : smear the <i>eT</i> according to a Gaussian
438with width <i>resolution * sqrt(eT)</i>, with the Gaussian truncated
439at 0 and <i>upperCut * eT</i>, 
440<br/><code>argumentoption </code><strong> 2</strong> : smear the <i>e = eT * cosh(eta)</i> according
441to a Gaussian with width <i>resolution * sqrt(e)</i>, with the
442Gaussian truncated at 0 and <i>upperCut * e</i>. 
443 
444<br/><code>argument</code><strong> resolution </strong> (<code>default = <strong>0.5</strong></code>) :
445see above.
446 
447<br/><code>argument</code><strong> upperCut </strong> (<code>default = <strong>2.</strong></code>) :
448see above.
449 
450<br/><code>argument</code><strong> threshold </strong> (<code>default = <strong>0 GeV</strong></code>) :
451completely neglect all bins with an <i>eT &lt; threshold</i>.
452 
453<br/><code>argument</code><strong> rndmPtr </strong> (<code>default = <strong>0</strong></code>) :
454the random-number generator used to select the smearing described
455above. Must be handed in for smearing to be possible. If your
456<code>Pythia</code> class instance is named <code>pythia</code>,
457then <code>&pythia.rndm</code> would be the logical choice.
458 
459 
460
461<a name="method26"></a>
462<p/><strong>bool CellJet::analyze( const Event& event, double eTjetMin = 20., double coneRadius = 0.7, double eTseed = 1.5, ostream& os = cout) &nbsp;</strong> <br/>
463performs a jet finding analysis, where
464<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
465most likely the <code>pythia.event</code> one.
466 
467<br/><code>argument</code><strong> eTjetMin </strong> (<code>default = <strong>20. GeV</strong></code>) : 
468is the minimum transverse energy inside a cone for this to be
469accepted as a jet.
470 
471<br/><code>argument</code><strong> coneRadius </strong> (<code>default = <strong>0.7</strong></code>) : 
472 is the size of the cone in <i>(eta, phi)</i> space drawn around
473the geometric center of the jet.
474 
475<br/><code>argument</code><strong> eTseed </strong> (<code>default = <strong>1.5 GeV</strong></code>) : 
476the mimimum <i>eT</i> in a cell for this to be acceptable as
477the trial center of a jet.
478 
479<br/><code>argument</code><strong> os </strong> (<code>default = <strong>cout</strong></code>) :  is the output stream for
480error messages. (The method does not rely on the <code>Info</code>
481mchinery for error messages.)
482 
483<br/>If the routine returns <code>false</code> the analysis failed,
484but currently this is not foreseen ever to happen.
485 
486
487<p/>
488After the analysis has been performed, a few <code>CellJet</code> 
489class methods are available to return the result of the analysis:
490
491<a name="method27"></a>
492<p/><strong>int CellJet::size() &nbsp;</strong> <br/>
493gives the number of jets found, with jets numbered 0 through
494<code>size() - 1</code>,
495 
496
497<a name="method28"></a>
498<p/><strong>double CellJet::eT(i) &nbsp;</strong> <br/>
499gives the <i>eT</i> of the <i>i</i>'th jet, where jets have been
500ordered with decreasing <i>eT</i> values,
501 
502
503<a name="method29"></a>
504<p/><strong>double CellJet::etaCenter(int i) &nbsp;</strong> <br/>
505 
506<strong>double CellJet::phiCenter(int i) &nbsp;</strong> <br/>
507gives the <i>eta</i> and <i>phi</i> coordinates of the geometrical
508center of the <i>i</i>'th jet,
509 
510
511<a name="method30"></a>
512<p/><strong>double CellJet::etaWeighted(int i) &nbsp;</strong> <br/>
513 
514<strong>double CellJet::phiWeighted(int i) &nbsp;</strong> <br/>
515gives the <i>eta</i> and <i>phi</i> coordinates of the
516<i>eT</i>-weighted center of the <i>i</i>'th jet,
517 
518
519<a name="method31"></a>
520<p/><strong>int CellJet::multiplicity(int i) &nbsp;</strong> <br/>
521gives the number of particles clustered into the <i>i</i>'th jet,
522 
523
524<a name="method32"></a>
525<p/><strong>Vec4 CellJet::pMassless(int i) &nbsp;</strong> <br/>
526gives a <code>Vec4</code> corresponding to the four-momentum defined
527by the <i>eT</i> and the weighted center of the <i>i</i>'th jet,
528 
529
530<a name="method33"></a>
531<p/><strong>Vec4 CellJet::pMassive(int i) &nbsp;</strong> <br/>
532gives a <code>Vec4</code> corresponding to the four-momentum defined by
533the sum of all the contributing cells to the <i>i</i>'th jet, where
534each cell contributes a four-momentum as if all the <i>eT</i> is
535deposited in the center of the cell,
536 
537
538<a name="method34"></a>
539<p/><strong>double CellJet::m(int i) &nbsp;</strong> <br/>
540gives the invariant mass of the <i>i</i>'th jet, defined by the
541<code>pMassive</code> above,
542 
543
544<a name="method35"></a>
545<p/><strong>void CellJet::list() &nbsp;</strong> <br/>
546provides a listing of the above information (except <code>pMassless</code>,
547for reasons of space).
548 
549
550<p/>
551There is also one method that returns information accumulated for all
552the events analyzed so far.
553<a name="method36"></a>
554<p/><strong>int CellJet::nError() &nbsp;</strong> <br/>
555tells the number of times <code>analyze(...)</code> failed to analyze
556events, i.e. returned <code>false</code>.
557 
558
559<h3>SlowJet</h3>
560
561<code>SlowJet</code> is a simple program for doing jet finding according
562to either of the <i>kT</i>, anti-<i>kT</i>, and Cambridge/Aachen
563algorithms, in a cylindrical coordinate frame. The name is obviously
564an homage to the <code>FastJet</code> program [<a href="Bibliography.html" target="page">Cac06</a>].
565That package contains many more algorithms, with many more options,
566and, above all, is <i>much</i> faster. Therefore <code>SlowJet</code>
567is not so much intended for massive processing of data or Monte Carlo
568files as for simple first tests. Nevertheless, within the advertised
569capabilities of <code>SlowJet</code>, it has been checked to find
570identically the same jets as <code>FastJet</code>. The time consumption
571typically is around or below that to generate an LHC <i>pp</i> event
572in the first place, so is not prohibitive. But the time rises rapidly
573for large multiplicities, so obviously <code>SlowJet</code> can not
574be used for tricks like distributing a dense grid of pseudoparticles
575to be able to define jet areas, like <code>FastJet</code> can, and also
576not for events with much pileup or other noise.
577
578<p/> 
579The first step is to decide which particles should be included in the
580analysis, and with what four-momenta. The <code>SlowJet</code> constructor
581allows to pick a maximum pseudorapidity defined by the extent of the
582assumed detector, to pick among some standard options of which particles
583to analyze, and to allow for some standard mass assumptions, like that
584all charged particles have the pion mass. Obviously this is only a
585restricted set of possibilities.
586
587<p/> 
588Full flexibility can be obtained by deriving from the base class
589<code>SlowJetHook</code> to write your own <code>include</code> method.
590This will be presented with one final-state particle at a time, and
591should return <code>true</code> for those particles that should be
592analyzed. It is also possible to return modified four-momenta and masses,
593to take into account detector smearing effects or particle identity
594misassignments, but you must respect <i>E^2 - p^2 = m^2</i>.
595
596<p/> 
597Alternatively you can modify the event record itself, or a copy of it
598(if you want to keep the original intact). For instance, only final
599particles are considered in the analysis, i.e. particles with positive
600status code, so negative status code should then be assigned to those
601particles that you do not want to see analyzed. Again four-momenta and
602masses can be modified, subject to <i>E^2 - p^2 = m^2</i>.
603
604<p/> 
605The jet reconstructions is then based on sequential recombination with
606progressive removal, using the <i>E</i> recombination scheme. To be
607more specific, the algorithm works as follows.
608<ol>
609<li>Each particle to be analyzed defines an original cluster. It has a
610well-defined four-momentum and mass at input. From this information the
611triplet <i>(pT, y, phi)</i> is calculated, i.e. the transverse momentum,
612rapidity and azimuthal angle of the cluster.</li>
613<li>Define distance measures of all clusters <i>i</i> to the beam
614<br/><i>d_iB = pT_i^2p</i><br/>
615and of all pairs <i>(i,j)</i> relative to each other
616<br/><i>d_ij = min( pT_i^2p, pT_j^2p) DeltaR_ij^2 / R^2 </i><br/>
617where
618<br/><i>DeltaR_ij^2 = (y_i - y_j)^2 + (phi_i - phi_j)^2.</i><br/>   
619The jet algorithm is determined by the user-selected <i>p</i> value,
620where <i>p = -1</i> corresponds to the anti-<i>kT</i> one,
621<i>p = 0</i> to the Cambridge/Aachen one and <i>p = 1</i> to the
622<i>kT</i> one. Also <i>R</i> is chosen by the user, to give an
623approximate measure of the size of jets. However, note that jets need
624not have a circular shape in <i>(y, phi)</i> space, so <i>R</i> 
625can not straight off be interpreted as a jet radius.</li>
626<li>Find the smallest of all <i>d_iB</i> and <i>d_ij</i>.</li>
627<li>If this is a <i>d_iB</i> then cluster <i>i</i> is removed from
628the clusters list and instead added to the jets list.
629Optionally, a <i>pTjetMin</i> requirement is imposed, where only
630clusters with <i>pT > pTjetMin</i> are added to the jets list.
631If so, some of the analyzed particles will not be part of any final
632jet.</li>
633<li>If instead te smallest measure is a <i>d_ij</i> then the
634four-momenta of the <i>i</i> and <i>j</i> clusters are added
635to define a single new cluster. Convert this four-momentum to a new
636<i>(pT, y, phi)</i> triplet and update the list of <i>d_iB</i> 
637and <i>d_ij</i>.</li>
638<li>Return to step 3 until no clusters remain.</li>
639</ol>
640
641<p/> 
642To do jet finding analyses you first have to set up a <code>SlowJet</code>
643instance, where the arguments of the constructor specifies the details
644of the subsequent analyses. Thereafter you can feed in events to it,
645one at a time, and have them analyzed by the <code>analyze</code> method.
646Information on the resulting jets can be extracted by a few different methods.
647The minimal procedure only requires one call per event to do the analysis.
648We will begin by presenting it, and only afterwards some extensions.
649
650<a name="method37"></a>
651<p/><strong>SlowJet::SlowJet(double power, double R, double pTjetMin = 0.,double etaMax = 25., int select = 2, int massSet = 2, SlowJetHook* sjHookPtr = 0) &nbsp;</strong> <br/>
652create a <code>SlowJet</code> instance, where
653<br/><code>argument</code><strong> power </strong>  : 
654tells (half) the power of the transverse-momentum dependence of the
655distance measure,
656<br/><code>argumentoption </code><strong> -1</strong> : the anti-<i>kT</i> algorithm, 
657<br/><code>argumentoption </code><strong> 0</strong> : the Cambridge/Aachen algorithm, and 
658<br/><code>argumentoption </code><strong> 1</strong> : the <i>kT</i> algorithm. 
659 
660<br/><code>argument</code><strong> R </strong>  : 
661the <i>R</i> size parameter, which is crudely related to the radius of
662the jet cone in <i>(y, phi)</i> space around the center of the jet.
663 
664<br/><code>argument</code><strong> pTjetMin </strong> (<code>default = <strong>0.0 GeV</strong></code>) :
665the minimum transverse momentum required for a cluster
666to become a jet. By default all clusters become jets, and therefore
667all analyzed particles are assigned to a jet.
668For comparisons with perturbative QCD, however, it is only meaningful
669to consider jets with a significant <i>pT</i>.
670 
671<br/><code>argument</code><strong> etaMax </strong> (<code>default = <strong>25.</strong></code>) : 
672the maximum +-pseudorapidity that the detector is assumed to cover.
673If you pick a value above 20 there is assumed to be full coverage
674(obviously only meaningful for theoretical studies).
675 
676<br/><code>argument</code><strong> select </strong> (<code>default = <strong>2</strong></code>) : 
677tells which particles are analyzed,
678<br/><code>argumentoption </code><strong> 1</strong> : all final-state particles, 
679<br/><code>argumentoption </code><strong> 2</strong> : all observable final-state particles,
680i.e. excluding neutrinos and other particles without strong or
681electromagnetic interactions (the <code>isVisible()</code> particle
682method),
683and 
684<br/><code>argumentoption </code><strong> 3</strong> : only charged final-state particles. 
685 
686<br/><code>argument</code><strong> massSet </strong> (<code>default = <strong>2</strong></code>) : masses assumed for the particles
687used in the analysis
688<br/><code>argumentoption </code><strong> 0</strong> : all massless, 
689<br/><code>argumentoption </code><strong> 1</strong> : photons are massless while all others are
690assigned the <i>pi+-</i> mass, and
691 
692<br/><code>argumentoption </code><strong> 2</strong> : all given their correct masses. 
693 
694<br/><code>argument</code><strong> sjHookPtr </strong> (<code>default = <strong>0</strong></code>) : 
695gives the possibility to send in your own selection routine for which
696particles should be part of the analysis; see further below on the
697<code>SlowJetHook</code> class. If this pointer is sent in nonzero,
698<code>etaMax</code> and <code>massSet</code> are disregarded,
699and <code>select</code> only gives the basic selection, to which
700the user can add further requirements.
701 
702 
703
704<a name="method38"></a>
705<p/><strong>bool SlowJet::analyze( const Event& event) &nbsp;</strong> <br/>
706performs a jet finding analysis, where
707<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
708most likely the <code>pythia.event</code> one.
709 
710<br/>If the routine returns <code>false</code> the analysis failed,
711but currently this is not foreseen ever to happen.
712 
713
714<p/>
715After the analysis has been performed, a few <code>SlowJet</code> 
716class methods are available to return the result of the analysis:
717
718<a name="method39"></a>
719<p/><strong>int SlowJet::sizeOrig() &nbsp;</strong> <br/>
720gives the original number of particles (and thus clusters) that the
721analysis starts with.
722 
723
724<a name="method40"></a>
725<p/><strong>int SlowJet::sizeJet() &nbsp;</strong> <br/>
726gives the number of jets found, with jets numbered 0 through
727<code>sizeJet() - 1</code>, and ordered in terms of decreasing
728transverse momentum values w.r.t. the beam axis,
729 
730
731<a name="method41"></a>
732<p/><strong>double SlowJet::pT(i) &nbsp;</strong> <br/>
733gives the transverse momentum <i>pT</i> of the <i>i</i>'th jet,
734 
735
736<a name="method42"></a>
737<p/><strong>double SlowJet::y(int i) &nbsp;</strong> <br/>
738 
739<strong>double SlowJet::phi(int i) &nbsp;</strong> <br/>
740gives the rapidity <i>y</i> and azimuthal angle <i>phi</i> 
741of the center of the <i>i</i>'th jet (defined by the vector sum
742of constituent four-momenta),
743 
744
745<a name="method43"></a>
746<p/><strong>Vec4 SlowJet::p(int i) &nbsp;</strong> <br/>
747 
748<strong>double SlowJet::m(int i) &nbsp;</strong> <br/>
749gives a <code>Vec4</code> corresponding to the four-momentum
750sum of the particles assigned to the <i>i</i>'th jet, and
751the invariant mass of this four-vector,
752 
753
754<a name="method44"></a>
755<p/><strong>int SlowJet::multiplicity(int i) &nbsp;</strong> <br/>
756gives the number of particles clustered into the <i>i</i>'th jet,
757 
758
759<a name="method45"></a>
760<p/><strong>int SlowJet::jetAssignment(int i) &nbsp;</strong> <br/>
761gives the index of the jet that the particle <i>i</i> of the event
762record belongs to, or -1 if there is no jet containing particle
763<i>i</i>,
764 
765
766<a name="method46"></a>
767<p/><strong>void SlowJet::removeJet(int i) &nbsp;</strong> <br/>
768removes the <i>i</i>'th jet,
769 
770
771<a name="method47"></a>
772<p/><strong>void SlowJet::list() &nbsp;</strong> <br/>
773provides a listing of the above information.
774 
775
776<p/> 
777These are the basic methods. For more sophisticated usage
778it is possible to trace the clustering, one step at a time. If so, the
779<code>setup</code> method should be used to read in the event and find
780the initial smallest distance. Each subsequent <code>doStep</code>
781will then do one cluster joining and find the new smallest distance.
782You can therefore interrogate which clusters will be joined next
783before the joining actually is performed. Alternatively you can take
784several steps in one go, or take steps down to a predetermined number
785of jets plus remaining clusters.
786
787<a name="method48"></a>
788<p/><strong>bool SlowJet::setup( const Event& event) &nbsp;</strong> <br/>
789selects the particles to be analyzed, calculates initial distances,
790and finds the initial smallest distance.
791<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
792most likely the <code>pythia.event</code> one.
793 
794<br/>If the routine returns <code>false</code> the setup failed,
795but currently this is not foreseen ever to happen.
796 
797
798<a name="method49"></a>
799<p/><strong>bool SlowJet::doStep() &nbsp;</strong> <br/>
800do the next step of the clustering. This can either be that two
801clusters are joined to one, or that a cluster is promoted to a jet
802(which is discarded if its <i>pT</i> value is below
803<code>pTjetMin</code>).
804<br/>The routine will only return <code>false</code> if there are no
805clusters left.
806 
807
808<a name="method50"></a>
809<p/><strong>bool SlowJet::doNSteps(int nStep) &nbsp;</strong> <br/>
810calls the <code>doStep()</code> method <code>nStep</code> times,
811if possible. Will return <code>false</code> if the list of clusters
812is emptied before then. The stored jet information is still perfectly
813fine; it is only the number of steps that is wrong.
814 
815
816<a name="method51"></a>
817<p/><strong>bool SlowJet::stopAtN(int nStop) &nbsp;</strong> <br/>
818calls the <code>doStep()</code> method until a total of <code>nStop</code>
819jet and cluster objects remain. Will return <code>false</code> if this
820is not possible, specifically if the number of objects already is smaller
821than <code>nStop</code> when the method is called. The stored jet and
822cluster information is still perfectly fine; it only does not have the
823expected multiplicity.
824 
825
826<a name="method52"></a>
827<p/><strong>int SlowJet::sizeAll() &nbsp;</strong> <br/>
828gives the total current number of jets and clusters. The jets are
829numbered 0 through <code>sizeJet() - 1</code>, while the clusters
830are numbered <code>sizeJet()</code> through <code>sizeAll() - 1</code>.
831(Internally jets and clusters are represented by two separate arrays,
832but are here presented in one flat range.) Note that the jets are ordered
833in decreasing <i>pT</i>, while the clusters are not ordered.
834 
835
836<p/>
837With this extension, the methods <code>double pT(int i)</code>,
838<code>double y(int i)</code>, <code>double phi(int i)</code>,
839<code>Vec4 p(int i)</code>, <code>double m(int i)</code> and
840<code>int multiplicity(int i)</code> can be used as before.
841Furthermore, <code>list()</code> generalizes
842
843<a name="method53"></a>
844<p/><strong>void SlowJet::list(bool listAll = false, ostream& os = cout) &nbsp;</strong> <br/>
845provides a listing of the above information.
846<br/><code>argument</code><strong> listAll </strong>  :  lists both jets and clusters if <code>true</code>,
847else only jets.
848 
849 
850 
851<p/>
852Three further methods can be used to check what will happen next.
853
854<a name="method54"></a>
855<p/><strong>int SlowJet::iNext() &nbsp;</strong> <br/>
856 
857<strong>int SlowJet::jNext() &nbsp;</strong> <br/>
858 
859<strong>double SlowJet::dNext() &nbsp;</strong> <br/>
860if the next step is to join two clusters, then the methods give
861the <i>(i,j, d_ij)</i> values, if instead to promote
862a cluster to a jet then <i>(i, -1, d_iB)</i>.
863If no clusters remain then <i>(-1, -1, 0.)</i>. Note that
864the cluster numbers are offset as described above, i.e. they begin at
865<code>sizeJet()</code>, which of course easily could be subtracted off.
866Also note that the jet and cluster lists are (moderately) reshuffled
867in each new step.
868 
869 
870<p/>
871Finally, and separately, the <code>SlowJetHook</code> class can be used
872for a more smart selection of which particles to include in the analysis.
873For instance, isolated and/or high-<i>pT</i> muons, electrons and
874photons should presumably be identified separately at an early stage,
875and then not clustered to jets.
876 
877<p/>
878Technically, it works like with <a href="UserHooks.html" target="page">User Hooks</a>.
879That is, PYTHIA contains the base class. You write a derived class.
880In the main program you create an instance of this class, and hand it
881in to <code>SlowJet</code>; in this case already as part of the
882constructor.
883
884<p/> 
885The following methods should be defined in your derived class.
886
887<a name="method55"></a>
888<p/><strong>SlowJetHook::SlowJetHook() &nbsp;</strong> <br/>
889 
890<strong>virtual SlowJetHook::~SlowJetHook() &nbsp;</strong> <br/>
891the constructor and destructor need not do anything, and if so you
892need not write your own destructor.
893 
894
895<a name="method56"></a>
896<p/><strong>virtual bool SlowJetHook::include(int iSel, const Event& event, Vec4& pSel, double& mSel) &nbsp;</strong> <br/>
897is the main method that you will need to write. It will be called
898once for each final-state particle in an event, subject to the
899value of the <code>select</code> switch in the <code>SlowJet</code>
900constructor. The value <code>select = 2</code> may be convenient
901since then you do not need to remove e.g. neutrinos yourself, but
902use <code>select = 1</code> for full control. The method should then
903return <code>true</code> if you want to see particle included in the
904jet clustering, and <code>false</code> if not.
905<br/><code>argument</code><strong> iSel </strong>  : is the index in the event record of the
906currently studied particle.
907 
908<br/><code>argument</code><strong> event </strong>  : is an object of the <code>Event</code> class,
909most likely the <code>pythia.event</code> one, where the currently
910studied particle is found.
911 
912<br/><code>argument</code><strong> pSel </strong>  :  is at input the four-momentum of the
913currently studied particle. You can change the values, e.g. to take
914into account energy smearing in the detector, to define the initial
915cluster value, without corrupting the event record itself.
916 
917<br/><code>argument</code><strong> mSel </strong>  :  is at input the mass of the currently studied
918particle. You can change the value, e.g. to take into account
919particle misidentification, to define the initial cluster value,
920without corrupting the event record itself. Note that the changes of
921<code>pSel</code> and <code>mSel</code> must be coordinated such that
922<i>E^2 - p^2 = m^2</i> holds.
923 
924 
925
926<p/> 
927It is also possible to define further methods of your own.
928One such could e.g. be called directly in the main program before the
929<code>analyze</code> method is called, to identify and bookkeep
930some event properties you may not want to reanalyze for each
931individual particle.
932
933</body>
934</html>
935
936<!-- Copyright (C) 2012 Torbjorn Sjostrand -->
Note: See TracBrowser for help on using the repository browser.