PYTHIA8/pythia8130/include/Analysis.h

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