SiStripTools/src/EventShape.cc

0001 #include "DPGAnalysis/SiStripTools/interface/EventShape.h"
0002 #include <DataFormats/TrackReco/interface/Track.h>
0003 #include <TVector3.h>
0004 #include <vector>
0005 #include <TMatrixDSym.h>
0006 #include <TMatrixDSymEigen.h>
0007 #include <TVectorD.h>
0008
0009 using namespace edm;
0010 using namespace reco;
0011 using namespace std;
0012 using reco::TrackCollection;
0013
0014 EventShape::EventShape(reco::TrackCollection& tracks) : eigenvalues(3) {
0015   for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
0016     p.push_back(TVector3(itTrack->px(), itTrack->py(), itTrack->pz()));
0017   }
0018
0019   // first fill the momentum tensor
0020   TMatrixDSym MomentumTensor(3);
0021   for (std::vector<TVector3>::const_iterator momentum = p.begin(); momentum < p.end(); ++momentum) {
0022     for (unsigned int i = 0; i < 3; i++)
0023       for (unsigned int j = 0; j <= i; j++) {
0024         MomentumTensor[i][j] += momentum[i] * momentum[j];
0025       }
0026   }
0027   MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
0028   // find the eigen values
0029   TMatrixDSymEigen eigen(MomentumTensor);
0030   TVectorD eigenvals = eigen.GetEigenValues();
0031   eigenvalues[0] = eigenvals[0];
0032   eigenvalues[1] = eigenvals[1];
0033   eigenvalues[2] = eigenvals[2];
0034   sort(eigenvalues.begin(), eigenvalues.end());
0035 }
0036
0037 math::XYZTLorentzVectorF EventShape::thrust() const {
0038   math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0, 0, 0, 0);
0039   TVector3 qtbo;
0040   TVector3 zero(0., 0., 0.);
0041   float vnew = 0.;
0042   uint32_t Np = p.size();
0043
0044   // for more than 2 tracks
0045   if (Np > 2) {
0046     float vmax = 0.;
0047     TVector3 vn, vm, vc, vl;
0048     for (unsigned int i = 0; i < Np - 1; i++)
0049       for (unsigned int j = i + 1; j < Np; j++) {
0050         vc = p[i].Cross(p[j]);
0051         vl = zero;
0052         for (unsigned int k = 0; k < Np; k++)
0053           if ((k != i) && (k != j)) {
0054             if (p[k].Dot(vc) >= 0.)
0055               vl = vl + p[k];
0056             else
0057               vl = vl - p[k];
0058           }
0059         // make all four sign-combinations for i,j
0060         vn = vl + p[j] + p[i];
0061         vnew = vn.Mag2();
0062         if (vnew > vmax) {
0063           vmax = vnew;
0064           vm = vn;
0065         }
0066         vn = vl + p[j] - p[i];
0067         vnew = vn.Mag2();
0068         if (vnew > vmax) {
0069           vmax = vnew;
0070           vm = vn;
0071         }
0072         vn = vl - p[j] + p[i];
0073         vnew = vn.Mag2();
0074         if (vnew > vmax) {
0075           vmax = vnew;
0076           vm = vn;
0077         }
0078         vn = vl - p[j] - p[i];
0079         vnew = vn.Mag2();
0080         if (vnew > vmax) {
0081           vmax = vnew;
0082           vm = vn;
0083         }
0084       }
0085     // sum momenta of all particles and iterate
0086     for (int iter = 1; iter <= 4; iter++) {
0087       qtbo = zero;
0088       for (unsigned int i = 0; i < Np; i++)
0089         if (vm.Dot(p[i]) >= 0.)
0090           qtbo = qtbo + p[i];
0091         else
0092           qtbo = qtbo - p[i];
0093       vnew = qtbo.Mag2();
0094       if (vnew == vmax)
0095         break;
0096       vmax = vnew;
0097       vm = qtbo;
0098     }
0099   }  // of if Np > 2
0100   else if (Np == 2)
0101     if (p[0].Dot(p[1]) >= 0.)
0102       qtbo = p[0] + p[1];
0103     else
0104       qtbo = p[0] - p[1];
0105   else if (Np == 1)
0106     qtbo = p[0];
0107   else {
0108     qtbo = zero;
0109     return output;
0110   }
0111   // normalize thrust -division by total momentum-
0112   float vsum = 0.;
0113   for (unsigned int i = 0; i < Np; i++)
0114     vsum = vsum + p[i].Mag();
0115   vnew = qtbo.Mag();
0116   float v = vnew / vsum;
0117   float x = qtbo.X() / vnew;
0118   float y = qtbo.Y() / vnew;
0119   float z = qtbo.Z() / vnew;
0120   output.SetPxPyPzE(x, y, z, v);
0121   return output;
0122 }
0123
0124 math::XYZTLorentzVectorF EventShape::thrust(const reco::TrackCollection& tracks) {
0125   std::vector<TVector3> pp;
0126   uint32_t Np = tracks.size();
0127   math::XYZTLorentzVectorF output = math::XYZTLorentzVectorF(0, 0, 0, 0);
0128   for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
0129     pp.push_back(TVector3(itTrack->px(), itTrack->py(), itTrack->pz()));
0130   }
0131   TVector3 qtbo;
0132   TVector3 zero(0., 0., 0.);
0133   float vnew = 0.;
0134
0135   // for more than 2 tracks
0136   if (Np > 2) {
0137     float vmax = 0.;
0138     TVector3 vn, vm, vc, vl;
0139     for (unsigned int i = 0; i < Np - 1; i++)
0140       for (unsigned int j = i + 1; j < Np; j++) {
0141         vc = pp[i].Cross(pp[j]);
0142         vl = zero;
0143         for (unsigned int k = 0; k < Np; k++)
0144           if ((k != i) && (k != j)) {
0145             if (pp[k].Dot(vc) >= 0.)
0146               vl = vl + pp[k];
0147             else
0148               vl = vl - pp[k];
0149           }
0150         // make all four sign-combinations for i,j
0151         vn = vl + pp[j] + pp[i];
0152         vnew = vn.Mag2();
0153         if (vnew > vmax) {
0154           vmax = vnew;
0155           vm = vn;
0156         }
0157         vn = vl + pp[j] - pp[i];
0158         vnew = vn.Mag2();
0159         if (vnew > vmax) {
0160           vmax = vnew;
0161           vm = vn;
0162         }
0163         vn = vl - pp[j] + pp[i];
0164         vnew = vn.Mag2();
0165         if (vnew > vmax) {
0166           vmax = vnew;
0167           vm = vn;
0168         }
0169         vn = vl - pp[j] - pp[i];
0170         vnew = vn.Mag2();
0171         if (vnew > vmax) {
0172           vmax = vnew;
0173           vm = vn;
0174         }
0175       }
0176     // sum momenta of all particles and iterate
0177     for (int iter = 1; iter <= 4; iter++) {
0178       qtbo = zero;
0179       for (unsigned int i = 0; i < Np; i++)
0180         if (vm.Dot(pp[i]) >= 0.)
0181           qtbo = qtbo + pp[i];
0182         else
0183           qtbo = qtbo - pp[i];
0184       vnew = qtbo.Mag2();
0185       if (vnew == vmax)
0186         break;
0187       vmax = vnew;
0188       vm = qtbo;
0189     }
0190   }  // of if Np > 2
0191   else if (Np == 2)
0192     if (pp[0].Dot(pp[1]) >= 0.)
0193       qtbo = pp[0] + pp[1];
0194     else
0195       qtbo = pp[0] - pp[1];
0196   else if (Np == 1)
0197     qtbo = pp[0];
0198   else {
0199     qtbo = zero;
0200     return output;
0201   }
0202   // normalize thrust -division by total momentum-
0203   float vsum = 0.;
0204   for (unsigned int i = 0; i < Np; i++)
0205     vsum = vsum + pp[i].Mag();
0206   vnew = qtbo.Mag();
0207   float v = vnew / vsum;
0208   float x = qtbo.X() / vnew;
0209   float y = qtbo.Y() / vnew;
0210   float z = qtbo.Z() / vnew;
0211   output.SetPxPyPzE(x, y, z, v);
0212   return output;
0213 }
0214
0215 float EventShape::sphericity(const reco::TrackCollection& tracks) {
0216   // a critical check
0217   if (tracks.empty())
0218     return 0;
0219
0220   // first fill the momentum tensor
0221   TMatrixDSym MomentumTensor(3);
0222   for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
0223     std::vector<double> momentum(3);
0224     momentum[0] = itTrack->px();
0225     momentum[1] = itTrack->py();
0226     momentum[2] = itTrack->pz();
0227     for (unsigned int i = 0; i < 3; i++)
0228       for (unsigned int j = 0; j <= i; j++) {
0229         MomentumTensor[i][j] += momentum[i] * momentum[j];
0230       }
0231   }
0232   MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
0233   // find the eigen values
0234   TMatrixDSymEigen eigen(MomentumTensor);
0235   TVectorD eigenvals = eigen.GetEigenValues();
0236   vector<float> eigenvaluess(3);
0237   eigenvaluess[0] = eigenvals[0];
0238   eigenvaluess[1] = eigenvals[1];
0239   eigenvaluess[2] = eigenvals[2];
0240   sort(eigenvaluess.begin(), eigenvaluess.end());
0241   // compute spericity
0242   float sph = (1.5 * (1 - eigenvaluess[2]));
0243   return sph;
0244 }
0245
0246 float EventShape::aplanarity(const reco::TrackCollection& tracks) {
0247   // a critical check
0248   if (tracks.empty())
0249     return 0;
0250   // first fill the momentum tensor
0251   TMatrixDSym MomentumTensor(3);
0252   for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
0253     std::vector<double> momentum(3);
0254     momentum[0] = itTrack->px();
0255     momentum[1] = itTrack->py();
0256     momentum[2] = itTrack->pz();
0257     for (unsigned int i = 0; i < 3; i++)
0258       for (unsigned int j = 0; j <= i; j++) {
0259         MomentumTensor[i][j] += momentum[i] * momentum[j];
0260       }
0261   }
0262   MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
0263   // find the eigen values
0264   TMatrixDSymEigen eigen(MomentumTensor);
0265   TVectorD eigenvals = eigen.GetEigenValues();
0266   vector<float> eigenvaluess(3);
0267   eigenvaluess[0] = eigenvals[0];
0268   eigenvaluess[1] = eigenvals[1];
0269   eigenvaluess[2] = eigenvals[2];
0270   sort(eigenvaluess.begin(), eigenvaluess.end());
0271   // compute aplanarity
0272   return (1.5 * eigenvaluess[0]);
0273 }
0274
0275 float EventShape::planarity(const reco::TrackCollection& tracks) {
0276   // First a critical check
0277   if (tracks.empty())
0278     return 0;
0279   // first fill the momentum tensor
0280   TMatrixDSym MomentumTensor(3);
0281   for (reco::TrackCollection::const_iterator itTrack = tracks.begin(); itTrack < tracks.end(); ++itTrack) {
0282     std::vector<double> momentum(3);
0283     momentum[0] = itTrack->px();
0284     momentum[1] = itTrack->py();
0285     momentum[2] = itTrack->pz();
0286     for (unsigned int i = 0; i < 3; i++)
0287       for (unsigned int j = 0; j <= i; j++) {
0288         MomentumTensor[i][j] += momentum[i] * momentum[j];
0289       }
0290   }
0291   MomentumTensor *= 1 / (MomentumTensor[0][0] + MomentumTensor[1][1] + MomentumTensor[2][2]);
0292   // find the eigen values
0293   TMatrixDSymEigen eigen(MomentumTensor);
0294   TVectorD eigenvals = eigen.GetEigenValues();
0295   vector<float> eigenvaluess(3);
0296   eigenvaluess[0] = eigenvals[0];
0297   eigenvaluess[1] = eigenvals[1];
0298   eigenvaluess[2] = eigenvals[2];
0299   sort(eigenvaluess.begin(), eigenvaluess.end());
0300   // compute planarity
0301   return (eigenvaluess[0] / eigenvaluess[1]);
0302 }
0303
0304 float EventShape::sphericity() const {
0305   // compute sphericity
0306   return (1.5 * (1 - eigenvalues[2]));
0307 }
0308
0309 float EventShape::aplanarity() const {
0310   // compute aplanarity
0311   return (1.5 * eigenvalues[0]);
0312 }
0313
0314 float EventShape::planarity() const {
0315   // compute planarity
0316   return (eigenvalues[0] / eigenvalues[1]);
0317 }