]>
Commit | Line | Data |
---|---|---|
5ad4eb21 | 1 | // Analysis.h is a part of the PYTHIA event generator. |
2 | // Copyright (C) 2008 Torbjorn Sjostrand. | |
3 | // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details. | |
4 | // Please respect the MCnet Guidelines, see GUIDELINES for details. | |
5 | ||
6 | // Header file for the Sphericity, Thrust, ClusterJet and CellJet classes. | |
7 | // Sphericity: sphericity analysis of the event. | |
8 | // Thrust: thrust analysis of the event. | |
9 | // ClusterJet: clustering jet finder. | |
10 | // CellJet: calorimetric cone jet finder. | |
11 | ||
12 | #ifndef Pythia8_Analysis_H | |
13 | #define Pythia8_Analysis_H | |
14 | ||
15 | #include "Basics.h" | |
16 | #include "Event.h" | |
17 | #include "PythiaStdlib.h" | |
18 | ||
19 | namespace Pythia8 { | |
20 | ||
21 | //************************************************************************** | |
22 | ||
23 | // Sphericity class. | |
24 | // This class performs (optionally modified) sphericity analysis on an event. | |
25 | ||
26 | class Sphericity { | |
27 | ||
28 | public: | |
29 | ||
30 | // Constructor. | |
31 | Sphericity(double powerIn = 2., int selectIn = 2) : power(powerIn), | |
32 | select(selectIn), nFew(0), nBack(0) {powerInt = 0; | |
33 | if (abs(power - 1.) < 0.01) powerInt = 1; | |
34 | if (abs(power - 2.) < 0.01) powerInt = 2; | |
35 | powerMod = 0.5 * power - 1.;} | |
36 | ||
37 | // Analyze event. | |
38 | bool analyze(const Event& event, ostream& os = cout); | |
39 | ||
40 | // Return info on results of analysis. | |
41 | double sphericity() const {return 1.5 * (eVal2 + eVal3);} | |
42 | double aplanarity() const {return 1.5 * eVal3;} | |
43 | double eigenValue(int i) const {return (i < 2) ? eVal1 : | |
44 | ( (i < 3) ? eVal2 : eVal3 ) ;} | |
45 | Vec4 eventAxis(int i) const {return (i < 2) ? eVec1 : | |
46 | ( (i < 3) ? eVec2 : eVec3 ) ;} | |
47 | ||
48 | // Provide a listing of the info. | |
49 | void list(ostream& os = cout); | |
50 | ||
51 | // Tell how many events could not be analyzed. | |
52 | int nError() const {return nFew + nBack;} | |
53 | ||
54 | private: | |
55 | ||
56 | // Constants: could only be changed in the code itself. | |
57 | static const int NSTUDYMIN, TIMESTOPRINT; | |
58 | static const double P2MIN, EIGENVALUEMIN; | |
59 | ||
60 | // Properties of analysis. | |
61 | double power; | |
62 | int select, powerInt; | |
63 | double powerMod; | |
64 | ||
65 | // Outcome of analysis. | |
66 | double eVal1, eVal2, eVal3; | |
67 | Vec4 eVec1, eVec2, eVec3; | |
68 | ||
69 | // Error statistics; | |
70 | int nFew, nBack; | |
71 | ||
72 | }; | |
73 | ||
74 | //************************************************************************** | |
75 | ||
76 | // Thrust class. | |
77 | // This class performs thrust analysis on an event. | |
78 | ||
79 | class Thrust { | |
80 | ||
81 | public: | |
82 | ||
83 | // Constructor. | |
84 | Thrust(int selectIn = 2) : select(selectIn), nFew(0) {} | |
85 | ||
86 | // Analyze event. | |
87 | bool analyze(const Event& event, ostream& os = cout); | |
88 | ||
89 | // Return info on results of analysis. | |
90 | double thrust() const {return eVal1;} | |
91 | double tMajor() const {return eVal2;} | |
92 | double tMinor() const {return eVal3;} | |
93 | double oblateness() const {return eVal2 - eVal3;} | |
94 | Vec4 eventAxis(int i) const {return (i < 2) ? eVec1 : | |
95 | ( (i < 3) ? eVec2 : eVec3 ) ;} | |
96 | ||
97 | // Provide a listing of the info. | |
98 | void list(ostream& os = cout); | |
99 | ||
100 | // Tell how many events could not be analyzed. | |
101 | int nError() const {return nFew;} | |
102 | ||
103 | private: | |
104 | ||
105 | // Constants: could only be changed in the code itself. | |
106 | static const int NSTUDYMIN, TIMESTOPRINT; | |
107 | static const double MAJORMIN; | |
108 | ||
109 | // Properties of analysis. | |
110 | int select; | |
111 | ||
112 | // Outcome of analysis. | |
113 | double eVal1, eVal2, eVal3; | |
114 | Vec4 eVec1, eVec2, eVec3; | |
115 | ||
116 | // Error statistics; | |
117 | int nFew; | |
118 | ||
119 | }; | |
120 | ||
121 | //************************************************************************** | |
122 | ||
123 | // SingleClusterJet class. | |
124 | // Simple helper class to ClusterJet for a jet and its contents. | |
125 | ||
126 | class SingleClusterJet { | |
127 | ||
128 | public: | |
129 | ||
130 | // Constructors. | |
131 | SingleClusterJet(Vec4 pJetIn = 0., int motherIn = 0) : | |
132 | pJet(pJetIn), mother(motherIn), daughter(0), multiplicity(1), | |
133 | isAssigned(false) {pAbs = max( PABSMIN, pJet.pAbs());} | |
134 | SingleClusterJet& operator=(const SingleClusterJet& j) { if (this != &j) | |
135 | { pJet = j.pJet; mother = j.mother; daughter = j.daughter; | |
136 | multiplicity = j.multiplicity; pAbs = j.pAbs; | |
137 | isAssigned = j.isAssigned;} return *this; } | |
138 | ||
139 | // Properties of jet. | |
140 | // Note: mother, daughter and isAssigned only used for original | |
141 | // particles, multiplicity and pTemp only for reconstructed jets. | |
142 | Vec4 pJet; | |
143 | int mother, daughter, multiplicity; | |
144 | bool isAssigned; | |
145 | double pAbs; | |
146 | Vec4 pTemp; | |
147 | ||
148 | // Distance measures (Lund, JADE, Durham) with friend. | |
149 | friend double dist2Fun(int measure, const SingleClusterJet& j1, | |
150 | const SingleClusterJet& j2); | |
151 | ||
152 | private: | |
153 | ||
154 | // Constants: could only be changed in the code itself. | |
155 | static const double PABSMIN; | |
156 | ||
157 | } ; | |
158 | ||
159 | //************************************************************************** | |
160 | ||
161 | // ClusterJet class. | |
162 | // This class performs a jet clustering according to different | |
163 | // distance measures: Lund, JADE or Durham. | |
164 | ||
165 | class ClusterJet { | |
166 | ||
167 | public: | |
168 | ||
169 | // Constructor. | |
170 | ClusterJet(string measureIn = "Lund", int selectIn = 2, int massSetIn = 2, | |
171 | bool preclusterIn = false, bool reassignIn = false) : measure(1), | |
172 | select(selectIn), massSet(massSetIn), doPrecluster(preclusterIn), | |
173 | doReassign(reassignIn), nFew(0) { | |
174 | char firstChar = toupper(measureIn[0]); | |
175 | if (firstChar == 'J') measure = 2; | |
176 | if (firstChar == 'D') measure = 3; | |
177 | piMass = ParticleDataTable::m0(211); | |
178 | } | |
179 | ||
180 | // Analyze event. | |
181 | bool analyze(const Event& event, double yScaleIn, double pTscaleIn, | |
182 | int nJetMinIn = 1, int nJetMaxIn = 0, ostream& os = cout); | |
183 | ||
184 | // Return info on jets produced. | |
185 | int size() const {return jets.size();} | |
186 | Vec4 p(int j) const {return jets[j].pJet;} | |
187 | ||
188 | // Return belonging of particle to one of the jets (-1 if none). | |
189 | int jetAssignment(int i) const { | |
190 | for (int iP = 0; iP < int(particles.size()); ++iP) | |
191 | if (particles[iP].mother == i) return particles[iP].daughter; | |
192 | return -1;} | |
193 | ||
194 | // Provide a listing of the info. | |
195 | void list(ostream& os = cout); | |
196 | ||
197 | // Tell how many events could not be analyzed. | |
198 | int nError() const {return nFew;} | |
199 | ||
200 | private: | |
201 | ||
202 | // Constants: could only be changed in the code itself. | |
203 | static const int TIMESTOPRINT; | |
204 | static const double PABSMIN, PRECLUSTERFRAC, PRECLUSTERSTEP; | |
205 | ||
206 | // Properties of analysis. | |
207 | int measure, select, massSet; | |
208 | bool doPrecluster, doReassign; | |
209 | double yScale, pTscale; | |
210 | int nJetMin, nJetMax; | |
211 | ||
212 | // Temporary results. | |
213 | double piMass, dist2Join, dist2BigMin, distPre, dist2Pre; | |
214 | vector<SingleClusterJet> particles; | |
215 | int nParticles; | |
216 | ||
217 | // Error statistics; | |
218 | int nFew; | |
219 | ||
220 | // Member functions for some operations (for clarity). | |
221 | void precluster(); | |
222 | void reassign(); | |
223 | ||
224 | // Outcome of analysis: ET-ordered list of jets. | |
225 | vector<SingleClusterJet> jets; | |
226 | ||
227 | }; | |
228 | ||
229 | //************************************************************************** | |
230 | ||
231 | // SingleCell class. | |
232 | // Simple helper class to CellJet for a cell and its contents. | |
233 | ||
234 | class SingleCell { | |
235 | ||
236 | public: | |
237 | ||
238 | // Constructor. | |
239 | SingleCell(int iCellIn = 0, double etaCellIn = 0., double phiCellIn = 0., | |
240 | double eTcellIn = 0., int multiplicityIn = 0) : iCell(iCellIn), | |
241 | etaCell(etaCellIn), phiCell(phiCellIn), eTcell(eTcellIn), | |
242 | multiplicity(multiplicityIn), canBeSeed(true), isUsed(false), | |
243 | isAssigned(false) {} | |
244 | ||
245 | // Properties of cell. | |
246 | int iCell; | |
247 | double etaCell, phiCell, eTcell; | |
248 | int multiplicity; | |
249 | bool canBeSeed, isUsed, isAssigned; | |
250 | ||
251 | } ; | |
252 | ||
253 | //************************************************************************** | |
254 | ||
255 | // SingleCellJet class. | |
256 | // Simple helper class to CellJet for a jet and its contents. | |
257 | ||
258 | class SingleCellJet { | |
259 | ||
260 | public: | |
261 | ||
262 | // Constructor. | |
263 | SingleCellJet(double eTjetIn = 0., double etaCenterIn = 0., | |
264 | double phiCenterIn = 0., double etaWeightedIn = 0., | |
265 | double phiWeightedIn = 0., int multiplicityIn = 0, | |
266 | Vec4 pMassiveIn = 0.) : eTjet(eTjetIn), etaCenter(etaCenterIn), | |
267 | phiCenter(phiCenterIn), etaWeighted(etaWeightedIn), | |
268 | phiWeighted(phiWeightedIn), multiplicity(multiplicityIn), | |
269 | pMassive(pMassiveIn) {} | |
270 | ||
271 | // Properties of jet. | |
272 | double eTjet, etaCenter, phiCenter, etaWeighted, phiWeighted; | |
273 | int multiplicity; | |
274 | Vec4 pMassive; | |
275 | ||
276 | } ; | |
277 | ||
278 | //************************************************************************** | |
279 | ||
280 | // CellJet class. | |
281 | // This class performs a cone jet search in (eta, phi, E_T) space. | |
282 | ||
283 | class CellJet { | |
284 | ||
285 | public: | |
286 | ||
287 | // Constructor. | |
288 | CellJet(double etaMaxIn = 5., int nEtaIn = 50, int nPhiIn = 32, | |
289 | int selectIn = 2, int smearIn = 0, double resolutionIn = 0.5, | |
290 | double upperCutIn = 2., double thresholdIn = 0.) : etaMax(etaMaxIn), | |
291 | nEta(nEtaIn), nPhi(nPhiIn), select(selectIn), smear(smearIn), | |
292 | resolution(resolutionIn), upperCut(upperCutIn), | |
293 | threshold(thresholdIn), nFew(0) { } | |
294 | ||
295 | // Analyze event. | |
296 | bool analyze(const Event& event, double eTjetMinIn = 20., | |
297 | double coneRadiusIn = 0.7, double eTseedIn = 1.5, ostream& os = cout); | |
298 | ||
299 | // Return info on results of analysis. | |
300 | int size() const {return jets.size();} | |
301 | double eT(int i) const {return jets[i].eTjet;} | |
302 | double etaCenter(int i) const {return jets[i].etaCenter;} | |
303 | double phiCenter(int i) const {return jets[i].phiCenter;} | |
304 | double etaWeighted(int i) const {return jets[i].etaWeighted;} | |
305 | double phiWeighted(int i) const {return jets[i].phiWeighted;} | |
306 | double multiplicity(int i) const {return jets[i].multiplicity;} | |
307 | Vec4 pMassless(int i) const {return jets[i].eTjet * Vec4( | |
308 | cos(jets[i].phiWeighted), sin(jets[i].phiWeighted), | |
309 | sinh(jets[i].etaWeighted), cosh(jets[i].etaWeighted) );} | |
310 | Vec4 pMassive(int i) const {return jets[i].pMassive;} | |
311 | double m(int i) const {return jets[i].pMassive.mCalc();} | |
312 | ||
313 | // Provide a listing of the info. | |
314 | void list(ostream& os = cout); | |
315 | ||
316 | // Tell how many events could not be analyzed: so far never. | |
317 | int nError() const {return nFew;} | |
318 | ||
319 | private: | |
320 | ||
321 | // Constants: could only be changed in the code itself. | |
322 | static const int TIMESTOPRINT; | |
323 | ||
324 | // Properties of analysis. | |
325 | double etaMax; | |
326 | int nEta, nPhi, select, smear; | |
327 | double resolution, upperCut, threshold; | |
328 | double eTjetMin, coneRadius, eTseed; | |
329 | ||
330 | // Error statistics; | |
331 | int nFew; | |
332 | ||
333 | // Outcome of analysis: ET-ordered list of jets. | |
334 | vector<SingleCellJet> jets; | |
335 | ||
336 | }; | |
337 | ||
338 | //************************************************************************** | |
339 | ||
340 | } // end namespace Pythia8 | |
341 | ||
342 | #endif // end Pythia8_Analysis_H | |
343 |