PYTHIA8/pythia8145/include/Analysis.h

   1 // Analysis.h is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 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) const;
  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) const;
  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   }
 178
 179   // Analyze event.
 180   bool analyze(const Event& event, double yScaleIn, double pTscaleIn,
 181     int nJetMinIn = 1, int nJetMaxIn = 0, ostream& os = cout);
 182
 183   // Return info on jets produced.
 184   int    size() const {return jets.size();}
 185   Vec4 p(int i) const {return jets[i].pJet;}
 186
 187   // Return belonging of particle to one of the jets (-1 if none).
 188   int jetAssignment(int i) const {
 189     for (int iP = 0; iP < int(particles.size()); ++iP)
 190     if (particles[iP].mother == i) return particles[iP].daughter;
 191     return -1;}
 192
 193   // Provide a listing of the info.
 194   void list(ostream& os = cout) const;
 195
 196   // Return info on clustering values.
 197   int    distanceSize() const {return distances.size();}
 198   double distance(int i) const {
 199     return (i < distanceSize()) ? distances[i] : 0.; }
 200
 201   // Tell how many events could not be analyzed.
 202   int nError() const {return nFew;}
 203
 204 private:
 205
 206   // Constants: could only be changed in the code itself.
 207   static const int    TIMESTOPRINT;
 208   static const double PIMASS, PABSMIN, PRECLUSTERFRAC, PRECLUSTERSTEP;
 209
 210   // Properties of analysis.
 211   int    measure, select, massSet;
 212   bool   doPrecluster, doReassign;
 213   double yScale, pTscale;
 214   int    nJetMin, nJetMax;
 215
 216   // Temporary results.
 217   double piMass, dist2Join, dist2BigMin, distPre, dist2Pre;
 218   vector<SingleClusterJet> particles;
 219   int    nParticles;
 220
 221   // Error statistics;
 222   int    nFew;
 223
 224   // Member functions for some operations (for clarity).
 225   void precluster();
 226   void reassign();
 227
 228   // Outcome of analysis: ET-ordered list of jets.
 229   vector<SingleClusterJet> jets;
 230
 231   // Outcome of analysis: the distance values where the jets were merged.
 232   deque<double> distances;
 233
 234 };
 235
 236 //==========================================================================
 237
 238 // SingleCell class.
 239 // Simple helper class to CellJet for a cell and its contents.
 240
 241 class SingleCell {
 242
 243 public:
 244
 245   // Constructor.
 246   SingleCell(int iCellIn = 0, double etaCellIn = 0., double phiCellIn = 0.,
 247     double eTcellIn = 0., int multiplicityIn = 0) : iCell(iCellIn),
 248     etaCell(etaCellIn), phiCell(phiCellIn), eTcell(eTcellIn),
 249     multiplicity(multiplicityIn), canBeSeed(true), isUsed(false),
 250     isAssigned(false) {}
 251
 252   // Properties of cell.
 253   int    iCell;
 254   double etaCell, phiCell, eTcell;
 255   int    multiplicity;
 256   bool   canBeSeed, isUsed, isAssigned;
 257
 258 } ;
 259
 260 //==========================================================================
 261
 262 // SingleCellJet class.
 263 // Simple helper class to CellJet for a jet and its contents.
 264
 265 class SingleCellJet {
 266
 267 public:
 268
 269   // Constructor.
 270   SingleCellJet(double eTjetIn = 0., double etaCenterIn = 0.,
 271     double phiCenterIn = 0., double etaWeightedIn = 0.,
 272     double phiWeightedIn = 0., int multiplicityIn = 0,
 273     Vec4 pMassiveIn = 0.) : eTjet(eTjetIn), etaCenter(etaCenterIn),
 274     phiCenter(phiCenterIn), etaWeighted(etaWeightedIn),
 275     phiWeighted(phiWeightedIn), multiplicity(multiplicityIn),
 276     pMassive(pMassiveIn) {}
 277
 278   // Properties of jet.
 279   double eTjet, etaCenter, phiCenter, etaWeighted, phiWeighted;
 280   int    multiplicity;
 281   Vec4   pMassive;
 282
 283 } ;
 284
 285 //==========================================================================
 286
 287 // CellJet class.
 288 // This class performs a cone jet search in (eta, phi, E_T) space.
 289
 290 class CellJet {
 291
 292 public:
 293
 294   // Constructor.
 295   CellJet(double etaMaxIn = 5., int nEtaIn = 50, int nPhiIn = 32,
 296     int selectIn = 2, int smearIn = 0, double resolutionIn = 0.5,
 297     double upperCutIn = 2., double thresholdIn = 0., Rndm* rndmPtrIn = 0)
 298     : etaMax(etaMaxIn), nEta(nEtaIn), nPhi(nPhiIn), select(selectIn),
 299     smear(smearIn), resolution(resolutionIn), upperCut(upperCutIn),
 300     threshold(thresholdIn), nFew(0), rndmPtr(rndmPtrIn) { }
 301
 302   // Analyze event.
 303   bool analyze(const Event& event, double eTjetMinIn = 20.,
 304     double coneRadiusIn = 0.7, double eTseedIn = 1.5, ostream& os = cout);
 305
 306   // Return info on results of analysis.
 307   int    size()              const {return jets.size();}
 308   double eT(int i)           const {return jets[i].eTjet;}
 309   double etaCenter(int i)    const {return jets[i].etaCenter;}
 310   double phiCenter(int i)    const {return jets[i].phiCenter;}
 311   double etaWeighted(int i)  const {return jets[i].etaWeighted;}
 312   double phiWeighted(int i)  const {return jets[i].phiWeighted;}
 313   int    multiplicity(int i) const {return jets[i].multiplicity;}
 314   Vec4   pMassless(int i)    const {return jets[i].eTjet * Vec4(
 315            cos(jets[i].phiWeighted),  sin(jets[i].phiWeighted),
 316           sinh(jets[i].etaWeighted), cosh(jets[i].etaWeighted) );}
 317   Vec4   pMassive(int i)     const {return jets[i].pMassive;}
 318   double m(int i)            const {return jets[i].pMassive.mCalc();}
 319
 320   // Provide a listing of the info.
 321   void list(ostream& os = cout) const;
 322
 323   // Tell how many events could not be analyzed: so far never.
 324   int nError() const {return nFew;}
 325
 326 private:
 327
 328   // Constants: could only be changed in the code itself.
 329   static const int    TIMESTOPRINT;
 330
 331   // Properties of analysis.
 332   double etaMax;
 333   int    nEta, nPhi, select, smear;
 334   double resolution, upperCut, threshold;
 335   double eTjetMin, coneRadius, eTseed;
 336
 337   // Error statistics;
 338   int    nFew;
 339
 340   // Outcome of analysis: ET-ordered list of jets.
 341   vector<SingleCellJet> jets;
 342
 343   // Pointer to the random number generator (needed for energy smearing).
 344   Rndm* rndmPtr;
 345
 346 };
 347
 348 //==========================================================================
 349
 350 } // end namespace Pythia8
 351
 352 #endif // end Pythia8_Analysis_H
 353