StRoot  1
Analysis.h
1 // Analysis.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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 // SlowJet: recombination algorithm; lightweight version of FastJet.
12
13 #ifndef Pythia8_Analysis_H
14 #define Pythia8_Analysis_H
15
16 #include "Basics.h"
17 #include "Event.h"
18 #include "PythiaStdlib.h"
19
20 namespace Pythia8 {
21
22 //==========================================================================
23
24 // Sphericity class.
25 // This class performs (optionally modified) sphericity analysis on an event.
26
27 class Sphericity {
28
29 public:
30
31  // Constructor.
32  Sphericity(double powerIn = 2., int selectIn = 2) : power(powerIn),
33  select(selectIn), nFew(0), nBack(0) {powerInt = 0;
34  if (abs(power - 1.) < 0.01) powerInt = 1;
35  if (abs(power - 2.) < 0.01) powerInt = 2;
36  powerMod = 0.5 * power - 1.;}
37
38  // Analyze event.
39  bool analyze(const Event& event, ostream& os = cout);
40
41  // Return info on results of analysis.
42  double sphericity() const {return 1.5 * (eVal2 + eVal3);}
43  double aplanarity() const {return 1.5 * eVal3;}
44  double eigenValue(int i) const {return (i < 2) ? eVal1 :
45  ( (i < 3) ? eVal2 : eVal3 ) ;}
46  Vec4 eventAxis(int i) const {return (i < 2) ? eVec1 :
47  ( (i < 3) ? eVec2 : eVec3 ) ;}
48
49  // Provide a listing of the info.
50  void list(ostream& os = cout) const;
51
52  // Tell how many events could not be analyzed.
53  int nError() const {return nFew + nBack;}
54
55 private:
56
57  // Constants: could only be changed in the code itself.
58  static const int NSTUDYMIN, TIMESTOPRINT;
59  static const double P2MIN, EIGENVALUEMIN;
60
61  // Properties of analysis.
62  double power;
63  int select, powerInt;
64  double powerMod;
65
66  // Outcome of analysis.
67  double eVal1, eVal2, eVal3;
68  Vec4 eVec1, eVec2, eVec3;
69
70  // Error statistics;
71  int nFew, nBack;
72
73 };
74
75 //==========================================================================
76
77 // Thrust class.
78 // This class performs thrust analysis on an event.
79
80 class Thrust {
81
82 public:
83
84  // Constructor.
85  Thrust(int selectIn = 2) : select(selectIn), nFew(0) {}
86
87  // Analyze event.
88  bool analyze(const Event& event, ostream& os = cout);
89
90  // Return info on results of analysis.
91  double thrust() const {return eVal1;}
92  double tMajor() const {return eVal2;}
93  double tMinor() const {return eVal3;}
94  double oblateness() const {return eVal2 - eVal3;}
95  Vec4 eventAxis(int i) const {return (i < 2) ? eVec1 :
96  ( (i < 3) ? eVec2 : eVec3 ) ;}
97
98  // Provide a listing of the info.
99  void list(ostream& os = cout) const;
100
101  // Tell how many events could not be analyzed.
102  int nError() const {return nFew;}
103
104 private:
105
106  // Constants: could only be changed in the code itself.
107  static const int NSTUDYMIN, TIMESTOPRINT;
108  static const double MAJORMIN;
109
110  // Properties of analysis.
111  int select;
112
113  // Outcome of analysis.
114  double eVal1, eVal2, eVal3;
115  Vec4 eVec1, eVec2, eVec3;
116
117  // Error statistics;
118  int nFew;
119
120 };
121
122 //==========================================================================
123
124 // SingleClusterJet class.
125 // Simple helper class to ClusterJet for a jet and its contents.
126
128
129 public:
130
131  // Constructors.
132  SingleClusterJet(Vec4 pJetIn = 0., int motherIn = 0) :
133  pJet(pJetIn), mother(motherIn), daughter(0), multiplicity(1),
134  isAssigned(false) {pAbs = max( PABSMIN, pJet.pAbs());}
135  SingleClusterJet& operator=(const SingleClusterJet& j) { if (this != &j)
136  { pJet = j.pJet; mother = j.mother; daughter = j.daughter;
137  multiplicity = j.multiplicity; pAbs = j.pAbs;
138  isAssigned = j.isAssigned;} return *this; }
139
140  // Properties of jet.
141  // Note: mother, daughter and isAssigned only used for original
142  // particles, multiplicity and pTemp only for reconstructed jets.
143  Vec4 pJet;
144  int mother, daughter, multiplicity;
145  bool isAssigned;
146  double pAbs;
147  Vec4 pTemp;
148
149  // Distance measures (Lund, JADE, Durham) with friend.
150  friend double dist2Fun(int measure, const SingleClusterJet& j1,
151  const SingleClusterJet& j2);
152
153 private:
154
155  // Constants: could only be changed in the code itself.
156  static const double PABSMIN;
157
158 } ;
159
160 //--------------------------------------------------------------------------
161
162 // Namespace function declarations; friend of SingleClusterJet.
163
164 // Distance measures (Lund, JADE, Durham) with friend.
165 double dist2Fun(int measure, const SingleClusterJet& j1,
166  const SingleClusterJet& j2);
167
168 //==========================================================================
169
170 // ClusterJet class.
171 // This class performs a jet clustering according to different
172 // distance measures: Lund, JADE or Durham.
173
174 class ClusterJet {
175
176 public:
177
178  // Constructor.
179  ClusterJet(string measureIn = "Lund", int selectIn = 2, int massSetIn = 2,
180  bool preclusterIn = false, bool reassignIn = false) : measure(1),
181  select(selectIn), massSet(massSetIn), doPrecluster(preclusterIn),
182  doReassign(reassignIn), nFew(0) {
183  char firstChar = toupper(measureIn[0]);
184  if (firstChar == 'J') measure = 2;
185  if (firstChar == 'D') measure = 3;
186  }
187
188  // Analyze event.
189  bool analyze(const Event& event, double yScaleIn, double pTscaleIn,
190  int nJetMinIn = 1, int nJetMaxIn = 0, ostream& os = cout);
191
192  // Return info on jets produced.
193  int size() const {return jets.size();}
194  Vec4 p(int i) const {return jets[i].pJet;}
195  int mult(int i) const {return jets[i].multiplicity;}
196
197  // Return belonging of particle to one of the jets (-1 if none).
198  int jetAssignment(int i) const {
199  for (int iP = 0; iP < int(particles.size()); ++iP)
200  if (particles[iP].mother == i) return particles[iP].daughter;
201  return -1;}
202
203  // Provide a listing of the info.
204  void list(ostream& os = cout) const;
205
206  // Return info on clustering values.
207  int distanceSize() const {return distances.size();}
208  double distance(int i) const {
209  return (i < distanceSize()) ? distances[i] : 0.; }
210
211  // Tell how many events could not be analyzed.
212  int nError() const {return nFew;}
213
214 private:
215
216  // Constants: could only be changed in the code itself.
217  static const int TIMESTOPRINT;
218  static const double PIMASS, PABSMIN, PRECLUSTERFRAC, PRECLUSTERSTEP;
219
220  // Properties of analysis.
221  int measure, select, massSet;
222  bool doPrecluster, doReassign;
223  double yScale, pTscale;
224  int nJetMin, nJetMax;
225
226  // Temporary results.
227  double piMass, dist2Join, dist2BigMin, distPre, dist2Pre;
228  vector<SingleClusterJet> particles;
229  int nParticles;
230
231  // Error statistics;
232  int nFew;
233
234  // Member functions for some operations (for clarity).
235  void precluster();
236  void reassign();
237
238  // Outcome of analysis: ET-ordered list of jets.
239  vector<SingleClusterJet> jets;
240
241  // Outcome of analysis: the distance values where the jets were merged.
242  deque<double> distances;
243
244 };
245
246 //==========================================================================
247
248 // SingleCell class.
249 // Simple helper class to CellJet for a cell and its contents.
250
251 class SingleCell {
252
253 public:
254
255  // Constructor.
256  SingleCell(int iCellIn = 0, double etaCellIn = 0., double phiCellIn = 0.,
257  double eTcellIn = 0., int multiplicityIn = 0) : iCell(iCellIn),
258  etaCell(etaCellIn), phiCell(phiCellIn), eTcell(eTcellIn),
259  multiplicity(multiplicityIn), canBeSeed(true), isUsed(false),
260  isAssigned(false) {}
261
262  // Properties of cell.
263  int iCell;
264  double etaCell, phiCell, eTcell;
265  int multiplicity;
266  bool canBeSeed, isUsed, isAssigned;
267
268 } ;
269
270 //==========================================================================
271
272 // SingleCellJet class.
273 // Simple helper class to CellJet for a jet and its contents.
274
276
277 public:
278
279  // Constructor.
280  SingleCellJet(double eTjetIn = 0., double etaCenterIn = 0.,
281  double phiCenterIn = 0., double etaWeightedIn = 0.,
282  double phiWeightedIn = 0., int multiplicityIn = 0,
283  Vec4 pMassiveIn = 0.) : eTjet(eTjetIn), etaCenter(etaCenterIn),
284  phiCenter(phiCenterIn), etaWeighted(etaWeightedIn),
285  phiWeighted(phiWeightedIn), multiplicity(multiplicityIn),
286  pMassive(pMassiveIn) {}
287
288  // Properties of jet.
289  double eTjet, etaCenter, phiCenter, etaWeighted, phiWeighted;
290  int multiplicity;
291  Vec4 pMassive;
292
293 } ;
294
295 //==========================================================================
296
297 // CellJet class.
298 // This class performs a cone jet search in (eta, phi, E_T) space.
299
300 class CellJet {
301
302 public:
303
304  // Constructor.
305  CellJet(double etaMaxIn = 5., int nEtaIn = 50, int nPhiIn = 32,
306  int selectIn = 2, int smearIn = 0, double resolutionIn = 0.5,
307  double upperCutIn = 2., double thresholdIn = 0., Rndm* rndmPtrIn = 0)
308  : etaMax(etaMaxIn), nEta(nEtaIn), nPhi(nPhiIn), select(selectIn),
309  smear(smearIn), resolution(resolutionIn), upperCut(upperCutIn),
310  threshold(thresholdIn), nFew(0), rndmPtr(rndmPtrIn) { }
311
312  // Analyze event.
313  bool analyze(const Event& event, double eTjetMinIn = 20.,
314  double coneRadiusIn = 0.7, double eTseedIn = 1.5, ostream& os = cout);
315
316  // Return info on results of analysis.
317  int size() const {return jets.size();}
318  double eT(int i) const {return jets[i].eTjet;}
319  double etaCenter(int i) const {return jets[i].etaCenter;}
320  double phiCenter(int i) const {return jets[i].phiCenter;}
321  double etaWeighted(int i) const {return jets[i].etaWeighted;}
322  double phiWeighted(int i) const {return jets[i].phiWeighted;}
323  int multiplicity(int i) const {return jets[i].multiplicity;}
324  Vec4 pMassless(int i) const {return jets[i].eTjet * Vec4(
325  cos(jets[i].phiWeighted), sin(jets[i].phiWeighted),
326  sinh(jets[i].etaWeighted), cosh(jets[i].etaWeighted) );}
327  Vec4 pMassive(int i) const {return jets[i].pMassive;}
328  double m(int i) const {return jets[i].pMassive.mCalc();}
329
330  // Provide a listing of the info.
331  void list(ostream& os = cout) const;
332
333  // Tell how many events could not be analyzed: so far never.
334  int nError() const {return nFew;}
335
336 private:
337
338  // Constants: could only be changed in the code itself.
339  static const int TIMESTOPRINT;
340
341  // Properties of analysis.
342  double etaMax;
343  int nEta, nPhi, select, smear;
344  double resolution, upperCut, threshold;
346
347  // Error statistics;
348  int nFew;
349
350  // Outcome of analysis: ET-ordered list of jets.
351  vector<SingleCellJet> jets;
352
353  // Pointer to the random number generator (needed for energy smearing).
354  Rndm* rndmPtr;
355
356 };
357
358 //==========================================================================
359
360 // SlowJetHook class.
361 // Base class, used to derive your own class with your selection criteria.
362
363 class SlowJetHook {
364
365 public:
366
367  // Destructor.
368  virtual ~SlowJetHook() { }
369
370  // Method to be overloaded.
371  // It will be called for all final-state particles, one at a time, and
372  // should return true if the particle should be analyzed, false if not.
373  // The particle is in location iSel of the event record.
374  // If you wish you can also modify the four-momentum and mass that will
375  // be used in the analysis, without affecting the event record itself,
376  // by changing pSel and mSel. Remember to respect E^2 - p^2 = m^2.
377  virtual bool include(int iSel, const Event& event, Vec4& pSel,
378  double& mSel) = 0;
379
380 };
381
382 //==========================================================================
383
384 // SingleSlowJet class.
385 // Simple helper class to SlowJet for a jet and its contents.
386
388
389 public:
390
391  // Constructors.
392  SingleSlowJet( Vec4 pIn = 0., double pT2In = 0., double yIn = 0.,
393  double phiIn = 0.) : p(pIn), pT2(pT2In), y(yIn), phi(phiIn),
394  mult(1) { }
395  SingleSlowJet(const SingleSlowJet& ssj) : p(ssj.p), pT2(ssj.pT2),
396  y(ssj.y), phi(ssj.phi), mult(ssj.mult) { }
397  SingleSlowJet& operator=(const SingleSlowJet& ssj) { if (this != &ssj)
398  { p = ssj.p; pT2 = ssj.pT2; y = ssj.y; phi = ssj.phi;
399  mult = ssj.mult;} return *this; }
400
401  // Properties of jet.
402  Vec4 p;
403  double pT2, y, phi;
404  int mult;
405
406 } ;
407
408 //==========================================================================
409
410 // SlowJet class.
411 // This class performs a recombination jet search in (y, phi, pT) space.
412
413 class SlowJet {
414
415 public:
416
417  // Constructor.
418  SlowJet(int powerIn, double Rin, double pTjetMinIn = 0.,
419  double etaMaxIn = 25., int selectIn = 2, int massSetIn = 2,
420  SlowJetHook* sjHookPtrIn = 0) : power(powerIn), R(Rin),
421  pTjetMin(pTjetMinIn), etaMax(etaMaxIn), select(selectIn),
422  massSet(massSetIn), sjHookPtr(sjHookPtrIn)
423  { isAnti = (power < 0); isKT = (power > 0); isCA = (power == 0);
424  R2 = R*R; pT2jetMin = pTjetMin*pTjetMin; cutInEta = (etaMax <= 20.);
425  chargedOnly = (select > 2); visibleOnly = (select == 2);
426  modifyMass = (massSet < 2); noHook = (sjHookPtr == 0); }
427
428  // Analyze event, all in one go.
429  bool analyze(const Event& event) {
430  if ( !setup(event) ) return false;
431  while (clSize > 0) doStep(); return true; }
432
433  // Set up list of particles to analyze, and initial distances.
434  bool setup(const Event& event);
435
436  // Do one recombination step, possibly giving a jet.
437  bool doStep();
438
439  // Do several recombinations steps, if possible.
440  bool doNSteps(int nStep) {
441  while(nStep > 0 && clSize > 0) { doStep(); --nStep;}
442  return (nStep == 0); }
443
444  // Do recombinations until fixed numbers of clusters and jets remain.
445  bool stopAtN(int nStop) {
446  while (clSize + jtSize > nStop && clSize > 0) doStep();
447  return (clSize + jtSize == nStop); }
448
449  // Return info on jet (+cluster) results of analysis.
450  int sizeOrig() const {return origSize;}
451  int sizeJet() const {return jtSize;}
452  int sizeAll() const {return jtSize + clSize;}
453  double pT(int i) const {return (i < jtSize)
454  ? sqrt(jets[i].pT2) : sqrt(clusters[i - jtSize].pT2);}
455  double y(int i) const {return (i < jtSize)
456  ? jets[i].y : clusters[i - jtSize].y;}
457  double phi(int i) const {return (i < jtSize)
458  ? jets[i].phi : clusters[i - jtSize].phi;}
459  Vec4 p(int i) const {return (i < jtSize)
460  ? jets[i].p : clusters[i - jtSize].p;}
461  double m(int i) const {return (i < jtSize)
462  ? jets[i].p.mCalc() : clusters[i - jtSize].p.mCalc();}
463  int multiplicity(int i) const {return (i < jtSize)
464  ? jets[i].mult : clusters[i - jtSize].mult;}
465
466  // Return info on next step to be taken.
467  int iNext() const {return (iMin == -1) ? -1 : iMin + jtSize;}
468  int jNext() const {return (jMin == -1) ? -1 : jMin + jtSize;}
469  double dNext() const {return dMin;}
470
471  // Provide a listing of the info.
472  void list(bool listAll = false, ostream& os = cout) const;
473
474 private:
475
476  // Constants: could only be changed in the code itself.
477  static const int TIMESTOPRINT;
478  static const double PIMASS, TINY;
479
480  // Properties of analysis as such.
481  int power;
482  double R, pTjetMin, etaMax, R2, pT2jetMin;
483  int select, massSet;
484  bool isAnti, isKT, isCA, cutInEta, chargedOnly, visibleOnly,
485  modifyMass, noHook;
486
487  // SlowJetHook can be use to tailor particle selection step.
488  SlowJetHook* sjHookPtr;
489
490  // Intermediate clustering objects and final jet objects.
491  vector<SingleSlowJet> clusters;
492  vector<SingleSlowJet> jets;
493
494  // Intermediate clustering distances.
495  vector<double> diB;
496  vector<double> dij;
497
498  // Other intermediate variables.
499  int origSize, clSize, clLast, jtSize, iMin, jMin;
500  double dPhi, dijTemp, dMin;
501
502  // Find next cluster pair to join.
503  void findNext();
504
505 };
506
507 //==========================================================================
508
509 } // end namespace Pythia8
510
511 #endif // end Pythia8_Analysis_H
512