CMS 3D CMS Logo

Public Member Functions | Static Public Member Functions | Private Attributes

EventShape Class Reference

#include <EventShape.h>

List of all members.

Public Member Functions

float aplanarity () const
 EventShape (reco::TrackCollection &)
float planarity () const
float sphericity () const
math::XYZTLorentzVectorF thrust () const

Static Public Member Functions

static float aplanarity (const reco::TrackCollection &)
static float planarity (const reco::TrackCollection &)
static float sphericity (const reco::TrackCollection &)
static math::XYZTLorentzVectorF thrust (const reco::TrackCollection &)

Private Attributes

std::vector< float > eigenvalues
std::vector< TVector3 > p

Detailed Description

Definition at line 6 of file EventShape.h.


Constructor & Destructor Documentation

EventShape::EventShape ( reco::TrackCollection tracks)

Definition at line 14 of file EventShape.cc.

References eigenvalues, i, j, p, and python::multivaluedict::sort().

                                                 :eigenvalues(3)
{
  for(reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack<tracks.end(); ++itTrack) {
    p.push_back(TVector3(itTrack->px(),itTrack->py(),itTrack->pz()));
  }
  
  // first fill the momentum tensor
  TMatrixDSym MomentumTensor(3);
  for(std::vector<TVector3>::const_iterator momentum = p.begin();momentum<p.end();++momentum) {
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<=i;j++) {
        MomentumTensor[i][j] += momentum[i]*momentum[j];
      }
  }
  MomentumTensor*=1/(MomentumTensor[0][0]+MomentumTensor[1][1]+MomentumTensor[2][2]);
  // find the eigen values
  TMatrixDSymEigen eigen(MomentumTensor);
  TVectorD eigenvals = eigen.GetEigenValues();
  eigenvalues[0] = eigenvals[0];
  eigenvalues[1] = eigenvals[1];
  eigenvalues[2] = eigenvals[2];
  sort(eigenvalues.begin(),eigenvalues.end());
}

Member Function Documentation

float EventShape::aplanarity ( const reco::TrackCollection tracks) [static]

Definition at line 243 of file EventShape.cc.

References i, j, and python::multivaluedict::sort().

Referenced by TrackerDpgAnalysis::analyze().

{
  // a critical check
  if (tracks.size()==0) return 0;
  // first fill the momentum tensor
  TMatrixDSym MomentumTensor(3);
  for(reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack<tracks.end(); ++itTrack) {
  std::vector<double> momentum(3);
  momentum[0] = itTrack->px();
  momentum[1] = itTrack->py();
  momentum[2] = itTrack->pz();
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<=i;j++) {
        MomentumTensor[i][j] += momentum[i]*momentum[j];
      }
  }
  MomentumTensor*=1/(MomentumTensor[0][0]+MomentumTensor[1][1]+MomentumTensor[2][2]);
  // find the eigen values
  TMatrixDSymEigen eigen(MomentumTensor);
  TVectorD eigenvals = eigen.GetEigenValues();
  vector<float> eigenvaluess(3);
  eigenvaluess[0] = eigenvals[0];
  eigenvaluess[1] = eigenvals[1];
  eigenvaluess[2] = eigenvals[2];
  sort(eigenvaluess.begin(),eigenvaluess.end());
  // compute aplanarity
  return ( 1.5*eigenvaluess[0]);
}
float EventShape::aplanarity ( ) const

Definition at line 307 of file EventShape.cc.

References eigenvalues.

{
  // compute aplanarity
  return ( 1.5*eigenvalues[0]);
}
float EventShape::planarity ( const reco::TrackCollection tracks) [static]

Definition at line 272 of file EventShape.cc.

References i, j, and python::multivaluedict::sort().

Referenced by TrackerDpgAnalysis::analyze().

{
  // First a critical check
  if (tracks.size()==0) return 0;
  // first fill the momentum tensor
  TMatrixDSym MomentumTensor(3);
  for(reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack<tracks.end(); ++itTrack) {
  std::vector<double> momentum(3);
  momentum[0] = itTrack->px();
  momentum[1] = itTrack->py();
  momentum[2] = itTrack->pz();
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<=i;j++) {
        MomentumTensor[i][j] += momentum[i]*momentum[j];
      }
  }
  MomentumTensor*=1/(MomentumTensor[0][0]+MomentumTensor[1][1]+MomentumTensor[2][2]);
  // find the eigen values
  TMatrixDSymEigen eigen(MomentumTensor);
  TVectorD eigenvals = eigen.GetEigenValues();
  vector<float> eigenvaluess(3);
  eigenvaluess[0] = eigenvals[0];
  eigenvaluess[1] = eigenvals[1];
  eigenvaluess[2] = eigenvals[2];
  sort(eigenvaluess.begin(),eigenvaluess.end());
  // compute planarity
  return (eigenvaluess[0]/eigenvaluess[1]);
}
float EventShape::planarity ( ) const

Definition at line 313 of file EventShape.cc.

References eigenvalues.

{
  // compute planarity
  return (eigenvalues[0]/eigenvalues[1]);
}
float EventShape::sphericity ( const reco::TrackCollection tracks) [static]

Definition at line 212 of file EventShape.cc.

References i, j, and python::multivaluedict::sort().

Referenced by TrackerDpgAnalysis::analyze().

{
  // a critical check
  if(tracks.size()==0) return 0;
  
  // first fill the momentum tensor
  TMatrixDSym MomentumTensor(3);
  for(reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack<tracks.end(); ++itTrack) {
  std::vector<double> momentum(3);
  momentum[0] = itTrack->px();
  momentum[1] = itTrack->py();
  momentum[2] = itTrack->pz();
    for(unsigned int i=0;i<3;i++)
      for(unsigned int j=0;j<=i;j++) {
        MomentumTensor[i][j] += momentum[i]*momentum[j];
      }
  }
  MomentumTensor*=1/(MomentumTensor[0][0]+MomentumTensor[1][1]+MomentumTensor[2][2]);
  // find the eigen values
  TMatrixDSymEigen eigen(MomentumTensor);
  TVectorD eigenvals = eigen.GetEigenValues();
  vector<float> eigenvaluess(3);
  eigenvaluess[0] = eigenvals[0];
  eigenvaluess[1] = eigenvals[1];
  eigenvaluess[2] = eigenvals[2];
  sort(eigenvaluess.begin(),eigenvaluess.end());
  // compute spericity
  float sph = ( 1.5*(1-eigenvaluess[2]));
  return sph;
}
float EventShape::sphericity ( ) const

Definition at line 301 of file EventShape.cc.

References eigenvalues.

{
  // compute sphericity
  return ( 1.5*(1-eigenvalues[2]));
}
math::XYZTLorentzVectorF EventShape::thrust ( const reco::TrackCollection tracks) [static]

Definition at line 123 of file EventShape.cc.

References i, j, gen::k, convertSQLitetoXML_cfg::output, createTree::pp, v, x, detailsBasic3DVector::y, z, and zero.

Referenced by TrackerDpgAnalysis::analyze().

{
  std::vector<TVector3> pp;   
  uint32_t Np = tracks.size();
  math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0,0,0,0);
  for(reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack<tracks.end(); ++itTrack) {
    pp.push_back(TVector3(itTrack->px(),itTrack->py(),itTrack->pz()));
  }
  TVector3 qtbo;
  TVector3 zero(0.,0.,0.);
  float vnew = 0.;
  
  // for more than 2 tracks
  if (Np > 2) { 
      float vmax  = 0.;
      TVector3 vn, vm, vc, vl;
      for(unsigned int i=0; i< Np-1; i++) 
        for(unsigned int j=i+1; j < Np; j++) {
            vc = pp[i].Cross(pp[j]);
            vl = zero; 
            for(unsigned int k=0; k<Np; k++)
              if ((k != i) && (k != j)) {
                if (pp[k].Dot(vc) >= 0.) vl = vl + pp[k];
                else vl = vl - pp[k];
              }
            // make all four sign-combinations for i,j
            vn = vl + pp[j] + pp[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl + pp[j] - pp[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl - pp[j] + pp[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl - pp[j] - pp[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }    
          }
      // sum momenta of all particles and iterate
      for(int iter=1; iter<=4; iter++) {  
          qtbo = zero;
          for(unsigned int i=0; i< Np; i++)
            if (vm.Dot(pp[i]) >= 0.) 
              qtbo = qtbo + pp[i];
            else 
              qtbo = qtbo - pp[i];
          vnew = qtbo.Mag2();
          if (vnew  == vmax) break;
          vmax = vnew;
          vm = qtbo;
        }
    }  // of if Np > 2
  else
    if (Np == 2)
      if (pp[0].Dot(pp[1]) >= 0.) 
        qtbo = pp[0] + pp[1];
      else
        qtbo = pp[0] - pp[1];
    else if (Np == 1)
      qtbo = pp[0];
    else {
        qtbo = zero;
        return output;
    }
  // normalize thrust -division by total momentum-
  float vsum = 0.;
  for(unsigned int i=0; i < Np; i++) vsum = vsum + pp[i].Mag();
  vnew  = qtbo.Mag();
  float v = vnew/vsum; 
  float x = qtbo.X()/vnew;
  float y = qtbo.Y()/vnew;
  float z = qtbo.Z()/vnew;
  output.SetPxPyPzE(x, y, z, v);
  return output;
}
math::XYZTLorentzVectorF EventShape::thrust ( ) const

Definition at line 38 of file EventShape.cc.

References i, j, gen::k, convertSQLitetoXML_cfg::output, p, v, x, detailsBasic3DVector::y, z, and zero.

{
  math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0,0,0,0);
  TVector3 qtbo;
  TVector3 zero(0.,0.,0.);
  float vnew = 0.;
  uint32_t Np = p.size();
  
  // for more than 2 tracks
  if (Np > 2) { 
      float vmax  = 0.;
      TVector3 vn, vm, vc, vl;
      for(unsigned int i=0; i< Np-1; i++) 
        for(unsigned int j=i+1; j < Np; j++) {
            vc = p[i].Cross(p[j]);
            vl = zero; 
            for(unsigned int k=0; k<Np; k++)
              if ((k != i) && (k != j)) {
                if (p[k].Dot(vc) >= 0.) vl = vl + p[k];
                else vl = vl - p[k];
              }
            // make all four sign-combinations for i,j
            vn = vl + p[j] + p[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl + p[j] - p[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl - p[j] + p[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }
            vn = vl - p[j] - p[i];
            vnew = vn.Mag2();
            if (vnew >  vmax) {  
                vmax = vnew;
                vm = vn;
            }    
          }
      // sum momenta of all particles and iterate
      for(int iter=1; iter<=4; iter++) {  
          qtbo = zero;
          for(unsigned int i=0; i< Np; i++)
            if (vm.Dot(p[i]) >= 0.) 
              qtbo = qtbo + p[i];
            else 
              qtbo = qtbo - p[i];
          vnew = qtbo.Mag2();
          if (vnew  == vmax) break;
          vmax = vnew;
          vm = qtbo;
        }
    }  // of if Np > 2
  else
    if (Np == 2)
      if (p[0].Dot(p[1]) >= 0.) 
        qtbo = p[0] + p[1];
      else
        qtbo = p[0] - p[1];
    else if (Np == 1)
      qtbo = p[0];
    else {
        qtbo = zero;
        return output;
    }
  // normalize thrust -division by total momentum-
  float vsum = 0.;
  for(unsigned int i=0; i < Np; i++) vsum = vsum + p[i].Mag();
  vnew  = qtbo.Mag();
  float v = vnew/vsum; 
  float x = qtbo.X()/vnew;
  float y = qtbo.Y()/vnew;
  float z = qtbo.Z()/vnew;
  output.SetPxPyPzE(x, y, z, v);
  return output;
}

Member Data Documentation

std::vector<float> EventShape::eigenvalues [private]

Definition at line 25 of file EventShape.h.

Referenced by aplanarity(), EventShape(), planarity(), and sphericity().

std::vector<TVector3> EventShape::p [private]

Definition at line 24 of file EventShape.h.

Referenced by EventShape(), and thrust().