ParticleFlowReco/src/PFDisplacedVertexCandidate.cc

0001 #include "DataFormats/ParticleFlowReco/interface/PFDisplacedVertexCandidate.h"
0002
0003 #include "DataFormats/TrackReco/interface/Track.h"
0004 #include "DataFormats/Math/interface/Point3D.h"
0005
0006 using namespace std;
0007 using namespace reco;
0008
0009 PFDisplacedVertexCandidate::PFDisplacedVertexCandidate() {}
0010
0011 void PFDisplacedVertexCandidate::addElement(const TrackBaseRef element) { elements_.push_back(element); }
0012
0013 void PFDisplacedVertexCandidate::setLink(
0014     const unsigned i1, const unsigned i2, const float dist, const GlobalPoint& dcaPoint, const VertexLinkTest test) {
0015   assert(test < LINKTEST_ALL);
0016
0017   unsigned index = 0;
0018   bool ok = matrix2vector(i1, i2, index);
0019
0020   if (ok) {
0021     //ignore the  -1, -1 pair
0022     if (dist > -0.5) {
0023       VertexLink& l = vertexLinkData_[index];
0024       l.distance_ = dist;
0025       l.dcaPoint_ = dcaPoint;
0026       l.test_ |= (1 << test);
0027     } else  //delete if existing
0028     {
0029       VertexLinkData::iterator it = vertexLinkData_.find(index);
0030       if (it != vertexLinkData_.end())
0031         vertexLinkData_.erase(it);
0032     }
0033
0034   } else {
0035     assert(0);
0036   }
0037 }
0038
0039 void PFDisplacedVertexCandidate::associatedElements(const unsigned i,
0040                                                     const VertexLinkData& vertexLinkData,
0041                                                     multimap<float, unsigned>& sortedAssociates,
0042                                                     const VertexLinkTest test) const {
0043   sortedAssociates.clear();
0044
0045   // i is too large
0046   if (i > elements_.size())
0047     return;
0048   // assert(i>=0); // i >= 0, since i is unsigned
0049   for (unsigned ie = 0; ie < elements_.size(); ie++) {
0050     // considered element itself
0051     if (ie == i) {
0052       continue;
0053     }
0054
0055     // Order the elements by increasing distance !
0056
0057     unsigned index = 0;
0058     if (!matrix2vector(i, ie, index))
0059       continue;
0060
0061     float c2 = -1;
0062     VertexLinkData::const_iterator it = vertexLinkData.find(index);
0063     if (it != vertexLinkData.end() && (((1 << test) & it->second.test_) != 0 || (test == LINKTEST_ALL)))
0064       c2 = it->second.distance_;
0065
0066     // not associated
0067     if (c2 < 0) {
0068       continue;
0069     }
0070
0071     sortedAssociates.insert(pair<float, unsigned>(c2, ie));
0072   }
0073 }
0074
0075 // -------- Provide useful information -------- //
0076
0077 PFDisplacedVertexCandidate::DistMap PFDisplacedVertexCandidate::r2Map() const {
0078   DistMap r2Map;
0079
0080   for (unsigned ie1 = 0; ie1 < elements_.size(); ie1++)
0081     for (unsigned ie2 = ie1 + 1; ie2 < elements_.size(); ie2++) {
0082       GlobalPoint P = dcaPoint(ie1, ie2);
0083       if (P.x() > 1e9)
0084         continue;
0085
0086       float r2 = P.x() * P.x() + P.y() * P.y() + P.z() * P.z();
0087
0088       r2Map.insert(pair<float, pair<int, int> >(r2, pair<int, int>(ie1, ie2)));
0089     }
0090
0091   return r2Map;
0092 }
0093
0094 PFDisplacedVertexCandidate::DistVector PFDisplacedVertexCandidate::r2Vector() const {
0095   DistVector r2Vector;
0096
0097   for (unsigned ie1 = 0; ie1 < elements_.size(); ie1++)
0098     for (unsigned ie2 = ie1 + 1; ie2 < elements_.size(); ie2++) {
0099       GlobalPoint P = dcaPoint(ie1, ie2);
0100       if (P.x() > 1e9)
0101         continue;
0102
0103       float r2 = P.x() * P.x() + P.y() * P.y() + P.z() * P.z();
0104
0105       r2Vector.push_back(r2);
0106     }
0107
0108   return r2Vector;
0109 }
0110
0111 PFDisplacedVertexCandidate::DistVector PFDisplacedVertexCandidate::distVector() const {
0112   DistVector distVector;
0113
0114   for (unsigned ie1 = 0; ie1 < elements_.size(); ie1++)
0115     for (unsigned ie2 = ie1 + 1; ie2 < elements_.size(); ie2++) {
0116       float d = dist(ie1, ie2);
0117       if (d < -0.5)
0118         continue;
0119
0120       distVector.push_back(d);
0121     }
0122
0123   return distVector;
0124 }
0125
0126 const GlobalPoint PFDisplacedVertexCandidate::dcaPoint(unsigned ie1, unsigned ie2) const {
0127   GlobalPoint dcaPoint(1e10, 1e10, 1e10);
0128
0129   unsigned index = 0;
0130   if (!matrix2vector(ie1, ie2, index))
0131     return dcaPoint;
0132   VertexLinkData::const_iterator it = vertexLinkData_.find(index);
0133   if (it != vertexLinkData_.end())
0134     dcaPoint = it->second.dcaPoint_;
0135
0136   return dcaPoint;
0137 }
0138
0139 // -------- Internal tools -------- //
0140
0141 bool PFDisplacedVertexCandidate::testLink(unsigned ie1, unsigned ie2) const {
0142   float d = dist(ie1, ie2);
0143   if (d < -0.5)
0144     return false;
0145   return true;
0146 }
0147
0148 const float PFDisplacedVertexCandidate::dist(unsigned ie1, unsigned ie2) const {
0149   float dist = -1;
0150
0151   unsigned index = 0;
0152   if (!matrix2vector(ie1, ie2, index))
0153     return dist;
0154   VertexLinkData::const_iterator it = vertexLinkData_.find(index);
0155   if (it != vertexLinkData_.end())
0156     dist = it->second.distance_;
0157
0158   return dist;
0159 }
0160
0161 // -------- Storage of the information -------- //
0162
0163 unsigned PFDisplacedVertexCandidate::vertexLinkDataSize() const {
0164   unsigned n = elements_.size();
0165
0166   // number of possible undirected links between n elements.
0167   // reflective links impossible.
0168
0169   return n * (n - 1) / 2;
0170 }
0171
0172 bool PFDisplacedVertexCandidate::matrix2vector(unsigned iindex, unsigned jindex, unsigned& index) const {
0173   unsigned size = elements_.size();
0174   if (iindex == jindex || iindex >= size || jindex >= size) {
0175     return false;
0176   }
0177
0178   if (iindex > jindex)
0179     swap(iindex, jindex);
0180
0181   index = jindex - iindex - 1;
0182
0183   if (iindex > 0) {
0184     index += iindex * size;
0185     unsigned missing = iindex * (iindex + 1) / 2;
0186     index -= missing;
0187   }
0188
0189   return true;
0190 }
0191
0192 void PFDisplacedVertexCandidate::Dump(ostream& out) const {
0193   if (!out)
0194     return;
0195
0196   const vector<TrackBaseRef>& elements = elements_;
0197   out << "\t--- DisplacedVertexCandidate ---  " << endl;
0198   out << "\tnumber of elements: " << elements.size() << endl;
0199
0200   // Build element label (string) : elid from type, layer and occurence number
0201   // use stringstream instead of sprintf to concatenate string and integer into string
0202   for (unsigned ie = 0; ie < elements.size(); ie++) {
0203     math::XYZPoint Pi(elements[ie].get()->innerPosition());
0204     math::XYZPoint Po(elements[ie].get()->outerPosition());
0205
0206     float innermost_radius = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y() + Pi.z() * Pi.z());
0207     float outermost_radius = sqrt(Po.x() * Po.x() + Po.y() * Po.y() + Po.z() * Po.z());
0208     float innermost_rho = sqrt(Pi.x() * Pi.x() + Pi.y() * Pi.y());
0209     float outermost_rho = sqrt(Po.x() * Po.x() + Po.y() * Po.y());
0210
0211     double pt = elements[ie]->pt();
0212
0213     out << "ie = " << elements[ie].key() << " pt = " << pt << " innermost hit radius = " << innermost_radius
0214         << " rho = " << innermost_rho << " outermost hit radius = " << outermost_radius << " rho = " << outermost_rho
0215         << endl;
0216   }
0217
0218   out << endl;
0219 }