VertexTools/src/SubsetHsmModeFinder3d.cc

0001 #include "RecoVertex/VertexTools/interface/hsm_3d.h"
0002 #include "RecoVertex/VertexTools/interface/SubsetHsmModeFinder3d.h"
0003
0004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0005
0006 #include <algorithm>
0007
0008 typedef std::pair<GlobalPoint, float> PointAndDistance;
0009
0010 namespace {
0011   struct compareByDistance {
0012     bool operator()(const PointAndDistance& p1, const PointAndDistance& p2) { return (p1.second < p2.second); };
0013   };
0014 }  // namespace
0015
0016 GlobalPoint SubsetHsmModeFinder3d::operator()(const std::vector<PointAndDistance>& values) const {
0017   if (values.empty()) {
0018     throw VertexException("SubsetHsmModeFinder3d: no value given.");
0019   };
0020
0021   std::vector<GlobalPoint> pts;
0022   pts.reserve(values.size() - 1);
0023   std::vector<PointAndDistance> sorted_values(values.size());
0024   partial_sort_copy(values.begin(), values.end(), sorted_values.begin(), sorted_values.end(), compareByDistance());
0025
0026   std::vector<PointAndDistance>::iterator end = sorted_values.end();
0027   std::vector<PointAndDistance>::iterator begin = sorted_values.begin();
0028
0029   float dmax = 0.004;  // 40 microns, as a first try.
0030
0031   // we want at least 30 values
0032   unsigned int min_num = values.size() < 30 ? values.size() : 30;
0033
0034   // we also want at least 50 % of all values
0035   if (values.size() > 2 * min_num)
0036     min_num = (int)values.size() / 2;
0037
0038   while (pts.size() < min_num) {
0039     // we cut at a dmax
0040     std::vector<PointAndDistance>::iterator i;
0041     for (i = begin; i != end && (i->second < dmax); ++i) {
0042       pts.push_back(i->first);
0043     };
0044     dmax += 0.003;  // add 30 microns with every iteration
0045     begin = i;
0046   };
0047
0048   GlobalPoint ret = hsm_3d(pts);
0049   return ret;
0050 }
0051
0052 SubsetHsmModeFinder3d* SubsetHsmModeFinder3d::clone() const { return new SubsetHsmModeFinder3d(*this); }