CMS 3D CMS Logo
List of all members | Public Member Functions | Static Public Member Functions | Private Attributes
EventShape Class Reference

#include <EventShape.h>

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::EventShape ( reco::TrackCollection &  tracks)

Definition at line 14 of file EventShape.cc.

14  : eigenvalues(3) {
15  for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
16  p.push_back(TVector3(itTrack->px(), itTrack->py(), itTrack->pz()));
17  }
18 
19  // first fill the momentum tensor
20  TMatrixDSym MomentumTensor(3);
21  for (std::vector<TVector3>::const_iterator momentum = p.begin(); momentum < p.end(); ++momentum) {
22  for (unsigned int i = 0; i < 3; i++)
23  for (unsigned int j = 0; j <= i; j++) {
24  MomentumTensor[i][j] += momentum[i] * momentum[j];
25  }
26  }
27  MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
28  // find the eigen values
29  TMatrixDSymEigen eigen(MomentumTensor);
30  TVectorD eigenvals = eigen.GetEigenValues();
31  eigenvalues[0] = eigenvals[0];
32  eigenvalues[1] = eigenvals[1];
33  eigenvalues[2] = eigenvals[2];
34  sort(eigenvalues.begin(), eigenvalues.end());
35 }

References eigenvalues, mps_fire::i, dqmiolumiharvest::j, p, jetUpdater_cfi::sort, and tracks.

Member Function Documentation

◆ aplanarity() [1/2]

float EventShape::aplanarity ( ) const

Definition at line 309 of file EventShape.cc.

309  {
310  // compute aplanarity
311  return (1.5 * eigenvalues[0]);
312 }

References eigenvalues.

◆ aplanarity() [2/2]

float EventShape::aplanarity ( const reco::TrackCollection &  tracks)
static

Definition at line 246 of file EventShape.cc.

246  {
247  // a critical check
248  if (tracks.empty())
249  return 0;
250  // first fill the momentum tensor
251  TMatrixDSym MomentumTensor(3);
252  for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
253  std::vector<double> momentum(3);
254  momentum[0] = itTrack->px();
255  momentum[1] = itTrack->py();
256  momentum[2] = itTrack->pz();
257  for (unsigned int i = 0; i < 3; i++)
258  for (unsigned int j = 0; j <= i; j++) {
259  MomentumTensor[i][j] += momentum[i] * momentum[j];
260  }
261  }
262  MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
263  // find the eigen values
264  TMatrixDSymEigen eigen(MomentumTensor);
265  TVectorD eigenvals = eigen.GetEigenValues();
266  vector<float> eigenvaluess(3);
267  eigenvaluess[0] = eigenvals[0];
268  eigenvaluess[1] = eigenvals[1];
269  eigenvaluess[2] = eigenvals[2];
270  sort(eigenvaluess.begin(), eigenvaluess.end());
271  // compute aplanarity
272  return (1.5 * eigenvaluess[0]);
273 }

References mps_fire::i, dqmiolumiharvest::j, jetUpdater_cfi::sort, and tracks.

Referenced by TrackerDpgAnalysis::analyze().

◆ planarity() [1/2]

float EventShape::planarity ( ) const

Definition at line 314 of file EventShape.cc.

314  {
315  // compute planarity
316  return (eigenvalues[0] / eigenvalues[1]);
317 }

References eigenvalues.

◆ planarity() [2/2]

float EventShape::planarity ( const reco::TrackCollection &  tracks)
static

Definition at line 275 of file EventShape.cc.

275  {
276  // First a critical check
277  if (tracks.empty())
278  return 0;
279  // first fill the momentum tensor
280  TMatrixDSym MomentumTensor(3);
281  for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
282  std::vector<double> momentum(3);
283  momentum[0] = itTrack->px();
284  momentum[1] = itTrack->py();
285  momentum[2] = itTrack->pz();
286  for (unsigned int i = 0; i < 3; i++)
287  for (unsigned int j = 0; j <= i; j++) {
288  MomentumTensor[i][j] += momentum[i] * momentum[j];
289  }
290  }
291  MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
292  // find the eigen values
293  TMatrixDSymEigen eigen(MomentumTensor);
294  TVectorD eigenvals = eigen.GetEigenValues();
295  vector<float> eigenvaluess(3);
296  eigenvaluess[0] = eigenvals[0];
297  eigenvaluess[1] = eigenvals[1];
298  eigenvaluess[2] = eigenvals[2];
299  sort(eigenvaluess.begin(), eigenvaluess.end());
300  // compute planarity
301  return (eigenvaluess[0] / eigenvaluess[1]);
302 }

References mps_fire::i, dqmiolumiharvest::j, jetUpdater_cfi::sort, and tracks.

Referenced by TrackerDpgAnalysis::analyze().

◆ sphericity() [1/2]

float EventShape::sphericity ( ) const

Definition at line 304 of file EventShape.cc.

304  {
305  // compute sphericity
306  return (1.5 * (1 - eigenvalues[2]));
307 }

References eigenvalues.

◆ sphericity() [2/2]

float EventShape::sphericity ( const reco::TrackCollection &  tracks)
static

Definition at line 215 of file EventShape.cc.

215  {
216  // a critical check
217  if (tracks.empty())
218  return 0;
219 
220  // first fill the momentum tensor
221  TMatrixDSym MomentumTensor(3);
222  for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
223  std::vector<double> momentum(3);
224  momentum[0] = itTrack->px();
225  momentum[1] = itTrack->py();
226  momentum[2] = itTrack->pz();
227  for (unsigned int i = 0; i < 3; i++)
228  for (unsigned int j = 0; j <= i; j++) {
229  MomentumTensor[i][j] += momentum[i] * momentum[j];
230  }
231  }
232  MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
233  // find the eigen values
234  TMatrixDSymEigen eigen(MomentumTensor);
235  TVectorD eigenvals = eigen.GetEigenValues();
236  vector<float> eigenvaluess(3);
237  eigenvaluess[0] = eigenvals[0];
238  eigenvaluess[1] = eigenvals[1];
239  eigenvaluess[2] = eigenvals[2];
240  sort(eigenvaluess.begin(), eigenvaluess.end());
241  // compute spericity
242  float sph = (1.5 * (1 - eigenvaluess[2]));
243  return sph;
244 }

References mps_fire::i, dqmiolumiharvest::j, jetUpdater_cfi::sort, and tracks.

Referenced by TrackerDpgAnalysis::analyze().

◆ thrust() [1/2]

math::XYZTLorentzVectorF EventShape::thrust ( ) const

Definition at line 37 of file EventShape.cc.

37  {
38  math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0, 0, 0, 0);
39  TVector3 qtbo;
40  TVector3 zero(0., 0., 0.);
41  float vnew = 0.;
42  uint32_t Np = p.size();
43 
44  // for more than 2 tracks
45  if (Np > 2) {
46  float vmax = 0.;
47  TVector3 vn, vm, vc, vl;
48  for (unsigned int i = 0; i < Np - 1; i++)
49  for (unsigned int j = i + 1; j < Np; j++) {
50  vc = p[i].Cross(p[j]);
51  vl = zero;
52  for (unsigned int k = 0; k < Np; k++)
53  if ((k != i) && (k != j)) {
54  if (p[k].Dot(vc) >= 0.)
55  vl = vl + p[k];
56  else
57  vl = vl - p[k];
58  }
59  // make all four sign-combinations for i,j
60  vn = vl + p[j] + p[i];
61  vnew = vn.Mag2();
62  if (vnew > vmax) {
63  vmax = vnew;
64  vm = vn;
65  }
66  vn = vl + p[j] - p[i];
67  vnew = vn.Mag2();
68  if (vnew > vmax) {
69  vmax = vnew;
70  vm = vn;
71  }
72  vn = vl - p[j] + p[i];
73  vnew = vn.Mag2();
74  if (vnew > vmax) {
75  vmax = vnew;
76  vm = vn;
77  }
78  vn = vl - p[j] - p[i];
79  vnew = vn.Mag2();
80  if (vnew > vmax) {
81  vmax = vnew;
82  vm = vn;
83  }
84  }
85  // sum momenta of all particles and iterate
86  for (int iter = 1; iter <= 4; iter++) {
87  qtbo = zero;
88  for (unsigned int i = 0; i < Np; i++)
89  if (vm.Dot(p[i]) >= 0.)
90  qtbo = qtbo + p[i];
91  else
92  qtbo = qtbo - p[i];
93  vnew = qtbo.Mag2();
94  if (vnew == vmax)
95  break;
96  vmax = vnew;
97  vm = qtbo;
98  }
99  } // of if Np > 2
100  else if (Np == 2)
101  if (p[0].Dot(p[1]) >= 0.)
102  qtbo = p[0] + p[1];
103  else
104  qtbo = p[0] - p[1];
105  else if (Np == 1)
106  qtbo = p[0];
107  else {
108  qtbo = zero;
109  return output;
110  }
111  // normalize thrust -division by total momentum-
112  float vsum = 0.;
113  for (unsigned int i = 0; i < Np; i++)
114  vsum = vsum + p[i].Mag();
115  vnew = qtbo.Mag();
116  float v = vnew / vsum;
117  float x = qtbo.X() / vnew;
118  float y = qtbo.Y() / vnew;
119  float z = qtbo.Z() / vnew;
120  output.SetPxPyPzE(x, y, z, v);
121  return output;
122 }

References mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, convertSQLitetoXML_cfg::output, p, findQualityFiles::v, x, y, z, and SiPixelPI::zero.

◆ thrust() [2/2]

math::XYZTLorentzVectorF EventShape::thrust ( const reco::TrackCollection &  tracks)
static

Definition at line 124 of file EventShape.cc.

124  {
125  std::vector<TVector3> pp;
126  uint32_t Np = tracks.size();
127  math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0, 0, 0, 0);
128  for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
129  pp.push_back(TVector3(itTrack->px(), itTrack->py(), itTrack->pz()));
130  }
131  TVector3 qtbo;
132  TVector3 zero(0., 0., 0.);
133  float vnew = 0.;
134 
135  // for more than 2 tracks
136  if (Np > 2) {
137  float vmax = 0.;
138  TVector3 vn, vm, vc, vl;
139  for (unsigned int i = 0; i < Np - 1; i++)
140  for (unsigned int j = i + 1; j < Np; j++) {
141  vc = pp[i].Cross(pp[j]);
142  vl = zero;
143  for (unsigned int k = 0; k < Np; k++)
144  if ((k != i) && (k != j)) {
145  if (pp[k].Dot(vc) >= 0.)
146  vl = vl + pp[k];
147  else
148  vl = vl - pp[k];
149  }
150  // make all four sign-combinations for i,j
151  vn = vl + pp[j] + pp[i];
152  vnew = vn.Mag2();
153  if (vnew > vmax) {
154  vmax = vnew;
155  vm = vn;
156  }
157  vn = vl + pp[j] - pp[i];
158  vnew = vn.Mag2();
159  if (vnew > vmax) {
160  vmax = vnew;
161  vm = vn;
162  }
163  vn = vl - pp[j] + pp[i];
164  vnew = vn.Mag2();
165  if (vnew > vmax) {
166  vmax = vnew;
167  vm = vn;
168  }
169  vn = vl - pp[j] - pp[i];
170  vnew = vn.Mag2();
171  if (vnew > vmax) {
172  vmax = vnew;
173  vm = vn;
174  }
175  }
176  // sum momenta of all particles and iterate
177  for (int iter = 1; iter <= 4; iter++) {
178  qtbo = zero;
179  for (unsigned int i = 0; i < Np; i++)
180  if (vm.Dot(pp[i]) >= 0.)
181  qtbo = qtbo + pp[i];
182  else
183  qtbo = qtbo - pp[i];
184  vnew = qtbo.Mag2();
185  if (vnew == vmax)
186  break;
187  vmax = vnew;
188  vm = qtbo;
189  }
190  } // of if Np > 2
191  else if (Np == 2)
192  if (pp[0].Dot(pp[1]) >= 0.)
193  qtbo = pp[0] + pp[1];
194  else
195  qtbo = pp[0] - pp[1];
196  else if (Np == 1)
197  qtbo = pp[0];
198  else {
199  qtbo = zero;
200  return output;
201  }
202  // normalize thrust -division by total momentum-
203  float vsum = 0.;
204  for (unsigned int i = 0; i < Np; i++)
205  vsum = vsum + pp[i].Mag();
206  vnew = qtbo.Mag();
207  float v = vnew / vsum;
208  float x = qtbo.X() / vnew;
209  float y = qtbo.Y() / vnew;
210  float z = qtbo.Z() / vnew;
211  output.SetPxPyPzE(x, y, z, v);
212  return output;
213 }

References mps_fire::i, dqmiolumiharvest::j, dqmdumpme::k, convertSQLitetoXML_cfg::output, createTree::pp, tracks, findQualityFiles::v, x, y, z, and SiPixelPI::zero.

Referenced by TrackerDpgAnalysis::analyze().

Member Data Documentation

◆ eigenvalues

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

Definition at line 22 of file EventShape.h.

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

◆ p

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

Definition at line 21 of file EventShape.h.

Referenced by EventShape(), Electron.Electron::ptErr(), and thrust().

DDAxes::y
mps_fire.i
i
Definition: mps_fire.py:428
convertSQLitetoXML_cfg.output
output
Definition: convertSQLitetoXML_cfg.py:72
math::XYZTLorentzVectorF
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< float > > XYZTLorentzVectorF
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:22
DDAxes::x
SiPixelPI::zero
Definition: SiPixelPayloadInspectorHelper.h:39
findQualityFiles.v
v
Definition: findQualityFiles.py:179
DDAxes::z
EventShape::p
std::vector< TVector3 > p
Definition: EventShape.h:21
dqmdumpme.k
k
Definition: dqmdumpme.py:60
tracks
const uint32_t *__restrict__ const HitContainer *__restrict__ TkSoA *__restrict__ tracks
Definition: CAHitNtupletGeneratorKernelsImpl.h:159
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
EventShape::eigenvalues
std::vector< float > eigenvalues
Definition: EventShape.h:22
createTree.pp
pp
Definition: createTree.py:17
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
Generated for CMSSW Reference Manual by doxygen 1.8.16