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