Line data Source code
1 : // Analysis.h is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "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, 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;}
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, 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 :
127 0 : class SingleClusterJet {
128 :
129 : public:
130 :
131 : // Constructors.
132 0 : SingleClusterJet(Vec4 pJetIn = 0., int motherIn = 0) :
133 0 : pJet(pJetIn), mother(motherIn), daughter(0), multiplicity(1),
134 0 : isAssigned(false) {pAbs = max( PABSMIN, pJet.pAbs());}
135 0 : SingleClusterJet& operator=(const SingleClusterJet& j) { if (this != &j)
136 0 : { pJet = j.pJet; mother = j.mother; daughter = j.daughter;
137 0 : multiplicity = j.multiplicity; pAbs = j.pAbs;
138 0 : 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 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 0 : double eTcellIn = 0., int multiplicityIn = 0) : iCell(iCellIn),
258 0 : etaCell(etaCellIn), phiCell(phiCellIn), eTcell(eTcellIn),
259 0 : multiplicity(multiplicityIn), canBeSeed(true), isUsed(false),
260 0 : 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 0 : 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 0 : Vec4 pMassiveIn = 0.) : eTjet(eTjetIn), etaCenter(etaCenterIn),
284 0 : phiCenter(phiCenterIn), etaWeighted(etaWeightedIn),
285 0 : phiWeighted(phiWeightedIn), multiplicity(multiplicityIn),
286 0 : 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;
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 0 : class SingleSlowJet {
388 :
389 : public:
390 :
391 : // Constructors.
392 0 : SingleSlowJet( Vec4 pIn = 0., double pT2In = 0., double yIn = 0.,
393 0 : double phiIn = 0., int idxIn = 0) : p(pIn), pT2(pT2In), y(yIn),
394 0 : phi(phiIn), mult(1) { idx.insert(idxIn); }
395 0 : SingleSlowJet(const SingleSlowJet& ssj) : p(ssj.p), pT2(ssj.pT2),
396 0 : y(ssj.y), phi(ssj.phi), mult(ssj.mult), idx(ssj.idx) { }
397 0 : SingleSlowJet& operator=(const SingleSlowJet& ssj) { if (this != &ssj)
398 0 : { p = ssj.p; pT2 = ssj.pT2; y = ssj.y; phi = ssj.phi;
399 0 : 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 : // Analyze event, all in one go.
431 : bool analyze(const Event& event) {
432 : if ( !setup(event) ) return false;
433 : if (useFJcore) return clusterFJ();
434 : while (clSize > 0) doStep(); return true; }
435 :
436 : // Set up list of particles to analyze, and initial distances.
437 : bool setup(const Event& event);
438 :
439 : // Do one recombination step, possibly giving a jet.
440 : bool doStep();
441 :
442 : // Do several recombinations steps, if possible.
443 : bool doNSteps(int nStep) { if (useFJcore) return false;
444 : while(nStep > 0 && clSize > 0) { doStep(); --nStep;}
445 : return (nStep == 0); }
446 :
447 : // Do recombinations until fixed numbers of clusters and jets remain.
448 : bool stopAtN(int nStop) { if (useFJcore) return false;
449 : while (clSize + jtSize > nStop && clSize > 0) doStep();
450 : return (clSize + jtSize == nStop); }
451 :
452 : // Return info on jet (+cluster) results of analysis.
453 : int sizeOrig() const {return origSize;}
454 : int sizeJet() const {return jtSize;}
455 : int sizeAll() const {return jtSize + clSize;}
456 : double pT(int i) const {return (i < jtSize)
457 : ? sqrt(jets[i].pT2) : sqrt(clusters[i - jtSize].pT2);}
458 : double y(int i) const {return (i < jtSize)
459 : ? jets[i].y : clusters[i - jtSize].y;}
460 : double phi(int i) const {return (i < jtSize)
461 : ? jets[i].phi : clusters[i - jtSize].phi;}
462 : Vec4 p(int i) const {return (i < jtSize)
463 : ? jets[i].p : clusters[i - jtSize].p;}
464 : double m(int i) const {return (i < jtSize)
465 : ? jets[i].p.mCalc() : clusters[i - jtSize].p.mCalc();}
466 : int multiplicity(int i) const {return (i < jtSize)
467 : ? jets[i].mult : clusters[i - jtSize].mult;}
468 :
469 : // Return info on next step to be taken.
470 : int iNext() const {return (iMin == -1) ? -1 : iMin + jtSize;}
471 : int jNext() const {return (jMin == -1) ? -1 : jMin + jtSize;}
472 : double dNext() const {return dMin;}
473 :
474 : // Provide a listing of the info.
475 : void list(bool listAll = false, ostream& os = cout) const;
476 :
477 : // Give a list of all particles in the jet.
478 : vector<int> constituents(int j) { vector<int> cTemp;
479 : for (set<int>::iterator idxTmp = jets[j].idx.begin();
480 : idxTmp != jets[j].idx.end(); ++idxTmp)
481 : cTemp.push_back( *idxTmp); return cTemp;
482 : }
483 :
484 : // Give a list of all particles in the cluster.
485 : vector<int> clusConstituents(int j) { vector<int> cTemp;
486 : for (set<int>::iterator idxTmp = clusters[j].idx.begin();
487 : idxTmp != clusters[j].idx.end(); ++idxTmp)
488 : cTemp.push_back( *idxTmp); return cTemp;
489 : }
490 :
491 : // Give the index of the jet that the particle i of the event record
492 : // belongs to. Returns -1 if particle i is not found in a jet.
493 : int jetAssignment(int i) {
494 : for (int j = 0; j < sizeJet(); ++j)
495 : if (jets[j].idx.find(i) != jets[j].idx.end()) return j;
496 : return -1;
497 : }
498 :
499 : // Remove a jet.
500 : void removeJet(int i) {
501 : if (i < 0 || i >= jtSize) return;
502 : jets.erase(jets.begin() + i);
503 : jtSize--;
504 : }
505 :
506 : private:
507 :
508 : // Constants: could only be changed in the code itself.
509 : static const int TIMESTOPRINT;
510 : static const double PIMASS, TINY;
511 :
512 : // Properties of analysis as such.
513 : int power;
514 : double R, pTjetMin, etaMax, R2, pT2jetMin;
515 : int select, massSet;
516 : // SlowJetHook can be used to tailor particle selection step.
517 : SlowJetHook* sjHookPtr;
518 : bool useFJcore, useStandardR, isAnti, isKT, cutInEta, chargedOnly,
519 : visibleOnly, modifyMass, noHook;
520 :
521 : // Intermediate clustering objects and final jet objects.
522 : vector<SingleSlowJet> clusters;
523 : vector<SingleSlowJet> jets;
524 :
525 : // Intermediate clustering distances.
526 : vector<double> diB;
527 : vector<double> dij;
528 :
529 : // Other intermediate variables.
530 : int origSize, clSize, clLast, jtSize, iMin, jMin;
531 : double dPhi, dijTemp, dMin;
532 :
533 : // Find next cluster pair to join.
534 : void findNext();
535 :
536 : // Use FJcore interface to perform clustering.
537 : bool clusterFJ();
538 :
539 : };
540 :
541 : //==========================================================================
542 :
543 : } // end namespace Pythia8
544 :
545 : #endif // end Pythia8_Analysis_H
|
