StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Analysis.h
1 // Analysis.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2018 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, 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 "Pythia8/Basics.h"
17 #include "Pythia8/Event.h"
18 #include "Pythia8/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) {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);
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() const;
51 
52  // Tell how many events could not be analyzed.
53  int nError() const {return nFew;}
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;
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);
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() 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 
127 class SingleClusterJet {
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);
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() 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 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 
275 class SingleCellJet {
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);
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() 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;
345  double eTjetMin, coneRadius, eTseed;
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 
387 class SingleSlowJet {
388 
389 public:
390 
391  // Constructors.
392  SingleSlowJet( Vec4 pIn = 0., double pT2In = 0., double yIn = 0.,
393  double phiIn = 0., int idxIn = 0) : p(pIn), pT2(pT2In), y(yIn),
394  phi(phiIn), mult(1) { idx.insert(idxIn); }
395  SingleSlowJet(const SingleSlowJet& ssj) : p(ssj.p), pT2(ssj.pT2),
396  y(ssj.y), phi(ssj.phi), mult(ssj.mult), idx(ssj.idx) { }
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; idx = ssj.idx; } return *this; }
400 
401  // Properties of jet.
402  Vec4 p;
403  double pT2, y, phi;
404  int mult;
405  set<int> idx;
406 
407 };
408 
409 //==========================================================================
410 
411 // SlowJet class.
412 // This class performs a recombination jet search in (y, phi, pT) space.
413 
414 class SlowJet {
415 
416 public:
417 
418  // Constructor.
419  SlowJet(int powerIn, double Rin, double pTjetMinIn = 0.,
420  double etaMaxIn = 25., int selectIn = 2, int massSetIn = 2,
421  SlowJetHook* sjHookPtrIn = 0, bool useFJcoreIn = true,
422  bool useStandardRin = true) : power(powerIn), R(Rin),
423  pTjetMin(pTjetMinIn), etaMax(etaMaxIn), select(selectIn),
424  massSet(massSetIn), sjHookPtr(sjHookPtrIn), useFJcore(useFJcoreIn),
425  useStandardR(useStandardRin) { isAnti = (power < 0); isKT = (power > 0);
426  R2 = R*R; pT2jetMin = pTjetMin*pTjetMin; cutInEta = (etaMax <= 20.);
427  chargedOnly = (select > 2); visibleOnly = (select == 2);
428  modifyMass = (massSet < 2); noHook = (sjHookPtr == 0);}
429 
430  // Destructor.
431  virtual ~SlowJet() {};
432 
433  // Analyze event, all in one go.
434  bool analyze(const Event& event) {
435  if ( !setup(event) ) return false;
436  if (useFJcore) return clusterFJ();
437  while (clSize > 0) doStep();
438  return true; }
439 
440  // Set up list of particles to analyze, and initial distances.
441  bool setup(const Event& event);
442 
443  // Do one recombination step, possibly giving a jet.
444  virtual bool doStep();
445 
446  // Do several recombinations steps, if possible.
447  bool doNSteps(int nStep) { if (useFJcore) return false;
448  while(nStep > 0 && clSize > 0) { doStep(); --nStep;}
449  return (nStep == 0); }
450 
451  // Do recombinations until fixed numbers of clusters and jets remain.
452  bool stopAtN(int nStop) { if (useFJcore) return false;
453  while (clSize + jtSize > nStop && clSize > 0) doStep();
454  return (clSize + jtSize == nStop); }
455 
456  // Return info on jet (+cluster) results of analysis.
457  int sizeOrig() const {return origSize;}
458  int sizeJet() const {return jtSize;}
459  int sizeAll() const {return jtSize + clSize;}
460  double pT(int i) const {return (i < jtSize)
461  ? sqrt(jets[i].pT2) : sqrt(clusters[i - jtSize].pT2);}
462  double y(int i) const {return (i < jtSize)
463  ? jets[i].y : clusters[i - jtSize].y;}
464  double phi(int i) const {return (i < jtSize)
465  ? jets[i].phi : clusters[i - jtSize].phi;}
466  Vec4 p(int i) const {return (i < jtSize)
467  ? jets[i].p : clusters[i - jtSize].p;}
468  double m(int i) const {return (i < jtSize)
469  ? jets[i].p.mCalc() : clusters[i - jtSize].p.mCalc();}
470  int multiplicity(int i) const {return (i < jtSize)
471  ? jets[i].mult : clusters[i - jtSize].mult;}
472 
473  // Return info on next step to be taken.
474  int iNext() const {return (iMin == -1) ? -1 : iMin + jtSize;}
475  int jNext() const {return (jMin == -1) ? -1 : jMin + jtSize;}
476  double dNext() const {return dMin;}
477 
478  // Provide a listing of the info.
479  void list(bool listAll = false) const;
480 
481  // Give a list of all particles in the jet.
482  vector<int> constituents(int j) { vector<int> cTemp;
483  for (set<int>::iterator idxTmp = jets[j].idx.begin();
484  idxTmp != jets[j].idx.end(); ++idxTmp) cTemp.push_back( *idxTmp);
485  return cTemp;
486  }
487 
488  // Give a list of all particles in the cluster.
489  vector<int> clusConstituents(int j) { vector<int> cTemp;
490  for (set<int>::iterator idxTmp = clusters[j].idx.begin();
491  idxTmp != clusters[j].idx.end(); ++idxTmp) cTemp.push_back( *idxTmp);
492  return cTemp;
493  }
494 
495  // Give the index of the jet that the particle i of the event record
496  // belongs to. Returns -1 if particle i is not found in a jet.
497  int jetAssignment(int i) {
498  for (int j = 0; j < sizeJet(); ++j)
499  if (jets[j].idx.find(i) != jets[j].idx.end()) return j;
500  return -1;
501  }
502 
503  // Remove a jet.
504  void removeJet(int i) {
505  if (i < 0 || i >= jtSize) return;
506  jets.erase(jets.begin() + i);
507  jtSize--;
508  }
509 
510 protected:
511 
512  // Constants: could only be changed in the code itself.
513  static const int TIMESTOPRINT;
514  static const double PIMASS, TINY;
515 
516  // Properties of analysis as such.
517  int power;
518  double R, pTjetMin, etaMax, R2, pT2jetMin;
519  int select, massSet;
520  // SlowJetHook can be used to tailor particle selection step.
521  SlowJetHook* sjHookPtr;
522  bool useFJcore, useStandardR, isAnti, isKT, cutInEta, chargedOnly,
523  visibleOnly, modifyMass, noHook;
524 
525  // Intermediate clustering objects and final jet objects.
526  vector<SingleSlowJet> clusters;
527  vector<SingleSlowJet> jets;
528 
529  // Intermediate clustering distances.
530  vector<double> diB;
531  vector<double> dij;
532 
533  // Other intermediate variables.
534  int origSize, clSize, clLast, jtSize, iMin, jMin;
535  double dPhi, dijTemp, dMin;
536 
537  // Find next cluster pair to join.
538  virtual void findNext();
539 
540  // Use FJcore interface to perform clustering.
541  bool clusterFJ();
542 
543 };
544 
545 //==========================================================================
546 
547 } // end namespace Pythia8
548 
549 #endif // end Pythia8_Analysis_H
Definition: AgUStep.h:26