CandUtils/interface/EventShapeVariables.h

0001 #ifndef EventShapeVariables_h
0002 #define EventShapeVariables_h
0003
0004 /**
0005    \class   EventShapeVariables EventShapeVariables.h "PhysicsTools/CandUtils/interface/EventShapeVariables.h"
0006
0007    \brief   Class for the calculation of several event shape variables
0008
0009    Class for the calculation of several event shape variables. Isotropy, sphericity,
0010    aplanarity and circularity are supported. The class supports vectors of 3d vectors
0011    and edm::Views of reco::Candidates as input. The 3d vectors can be given in
0012    cartesian, cylindrical or polar coordinates. It exploits the ROOT::TMatrixDSym
0013    for the calculation of the sphericity and aplanarity.
0014
0015    See https://arxiv.org/pdf/hep-ph/0603175v2.pdf#page=524
0016    for an explanation of sphericity, aplanarity and the quantities C and D.
0017
0018    Author: Sebastian Naumann-Emme, University of Hamburg
0019            Roger Wolf, University of Hamburg
0020        Christian Veelken, UC Davis
0021 */
0022
0023 #include "DataFormats/Math/interface/Vector3D.h"
0024 #include "DataFormats/Candidate/interface/Candidate.h"
0025 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0026
0027 #include "TMatrixDSym.h"
0028 #include "TVectorD.h"
0029
0030 #include <vector>
0031
0032 class EventShapeVariables {
0033 public:
0034   /// constructor from reco::Candidates
0035   explicit EventShapeVariables(const edm::View<reco::Candidate>& inputVectors);
0036   /// constructor from XYZ coordinates
0037   explicit EventShapeVariables(const std::vector<math::XYZVector>& inputVectors);
0038   /// constructor from rho eta phi coordinates
0039   explicit EventShapeVariables(const std::vector<math::RhoEtaPhiVector>& inputVectors);
0040   /// constructor from r theta phi coordinates
0041   explicit EventShapeVariables(const std::vector<math::RThetaPhiVector>& inputVectors);
0042   /// default destructor
0043   ~EventShapeVariables() {}
0044
0045   /// the return value is 1 for spherical events and 0 for events linear in r-phi. This function
0046   /// needs the number of steps to determine how fine the granularity of the algorithm in phi
0047   /// should be
0048   double isotropy(const unsigned int& numberOfSteps = 1000) const;
0049
0050   /// the return value is 1 for spherical and 0 linear events in r-phi. This function needs the
0051   /// number of steps to determine how fine the granularity of the algorithm in phi should be
0052   double circularity(const unsigned int& numberOfSteps = 1000) const;
0053
0054   /// set exponent for computation of momentum tensor and related products
0055   void set_r(double r);
0056
0057   /// set number of Fox-Wolfram moments to compute
0058   void setFWmax(unsigned m);
0059
0060   /// 1.5*(q1+q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor
0061   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return values are 1 for spherical, 3/4 for
0062   /// plane and 0 for linear events
0063   double sphericity();
0064   /// 1.5*q2 where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor
0065   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return values are 0.5 for spherical and 0
0066   /// for plane and linear events
0067   double aplanarity();
0068   /// 3.*(q0*q1+q0*q2+q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor
0069   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return value is between 0 and 1
0070   /// and measures the 3-jet structure of the event (C vanishes for a "perfect" 2-jet event)
0071   double C();
0072   /// 27.*(q0*q1*q2) where q0>=q1>=q2>=0 are the eigenvalues of the momentum tensor
0073   /// sum{p_j[a]*p_j[b]}/sum{p_j**2} normalized to 1. Return value is between 0 and 1
0074   /// and measures the 4-jet structure of the event (D vanishes for a planar event)
0075   double D();
0076
0077   const std::vector<double>& getEigenValues() {
0078     if (!tensors_computed_)
0079       compTensorsAndVectors();
0080     return eigenValues_;
0081   }
0082   const std::vector<double>& getEigenValuesNoNorm() {
0083     if (!tensors_computed_)
0084       compTensorsAndVectors();
0085     return eigenValuesNoNorm_;
0086   }
0087   const TMatrixD& getEigenVectors() {
0088     if (!tensors_computed_)
0089       compTensorsAndVectors();
0090     return eigenVectors_;
0091   }
0092
0093   double getFWmoment(unsigned l);
0094   const std::vector<double>& getFWmoments();
0095
0096 private:
0097   /// helper function to fill the 3 dimensional momentum tensor from the inputVectors where needed
0098   /// also fill the 3 dimensional vectors of eigen-values and eigen-vectors;
0099   /// the largest (smallest) eigen-value is stored at index position 0 (2)
0100   void compTensorsAndVectors();
0101
0102   void computeFWmoments();
0103
0104   /// caching of input vectors
0105   std::vector<math::XYZVector> inputVectors_;
0106
0107   /// caching of output
0108   double r_;
0109   bool tensors_computed_;
0110   TMatrixD eigenVectors_;
0111   TVectorD eigenValuesTmp_, eigenValuesNoNormTmp_;
0112   std::vector<double> eigenValues_, eigenValuesNoNorm_;
0113
0114   /// Owen ; save computed Fox-Wolfram moments
0115   unsigned fwmom_maxl_;
0116   std::vector<double> fwmom_;
0117   bool fwmom_computed_;
0118 };
0119
0120 #endif