VertexTools/src/SMS.cc

0001 #ifdef PROJECT_NAME
0002 #include "RecoVertex/VertexTools/interface/SMS.h"
0003 #else
0004 #include <vector>
0005 #include <algorithm>
0006 #include "RecoVertex/VertexTools/interface/SMS.h"
0007 using namespace std;
0008 #endif
0009
0010 namespace {
0011
0012   typedef std::pair<float, const GlobalPoint*> MyPair;
0013   typedef std::pair<float, float> FloatPair;
0014   typedef std::pair<GlobalPoint, float> GlPtWt;
0015   typedef std::pair<float, const GlPtWt*> MyPairWt;
0016
0017   struct Sorter {
0018     bool operator()(const MyPair& pair1, const MyPair& pair2) { return (pair1.first < pair2.first); };
0019
0020     bool operator()(const FloatPair& pair1, const FloatPair& pair2) { return (pair1.first < pair2.first); };
0021
0022     bool operator()(const MyPairWt& pair1, const MyPairWt& pair2) { return (pair1.first < pair2.first); };
0023   };
0024
0025   bool debug() { return false; }
0026
0027   inline GlobalPoint& operator+=(GlobalPoint& a, const GlobalPoint& b) {
0028     a = GlobalPoint(a.x() + b.x(), a.y() + b.y(), a.z() + b.z());
0029     return a;
0030   }
0031   inline GlobalPoint& operator/=(GlobalPoint& a, float b) {
0032     a = GlobalPoint(a.x() / b, a.y() / b, a.z() / b);
0033     return a;
0034   }
0035
0036   GlobalPoint average(const std::vector<MyPair>& pairs, int nq) {
0037     GlobalPoint location(0, 0, 0);
0038     for (std::vector<MyPair>::const_iterator i = pairs.begin(); i != (pairs.begin() + nq); ++i)
0039       location += *(i->second);
0040     location /= nq;
0041     return location;
0042   }
0043
0044   GlobalPoint average(const std::vector<MyPairWt>& pairs, int nq) {
0045     GlobalPoint location(0, 0, 0);
0046     for (std::vector<MyPairWt>::const_iterator i = pairs.begin(); i != (pairs.begin() + nq); ++i)
0047       location += (i->second)->first;
0048     location /= nq;
0049     return location;
0050   }
0051
0052   typedef SMS::SMSType SMSType;
0053 }  // namespace
0054
0055 SMS::SMS(SMSType tp, float q) : theType(tp), theRatio(q) {}
0056
0057 GlobalPoint SMS::location(const std::vector<GlobalPoint>& data) const {
0058   if (theType & Weighted) {
0059     std::cout << "[SMS] warning: Weighted SMS was asked for, but data are "
0060               << "weightless!" << std::endl;
0061   };
0062   int nobs = data.size();
0063   int nq = (int)ceil(theRatio * nobs);
0064   // cout << "nobs= " << nobs << "  nq= " << nq << endl;
0065
0066   // Compute distances
0067   std::vector<MyPair> pairs;
0068
0069   for (std::vector<GlobalPoint>::const_iterator i = data.begin(); i != data.end(); ++i) {
0070     std::vector<float> D;
0071     // Compute squared distances to all points
0072     for (std::vector<GlobalPoint>::const_iterator j = data.begin(); j != data.end(); ++j) {
0073       D.push_back((*j - *i).mag2());
0074     }
0075     // Find q-quantile in each row of the distance matrix
0076     sort(D.begin(), D.end());
0077     MyPair tmp(D[nq - 1], &(*i));
0078     pairs.push_back(tmp);
0079   };
0080
0081   // Sort pairs by first element
0082   sort(pairs.begin(), pairs.end(), Sorter());
0083   if (!(theType & SMS::Interpolate) && !(theType & SMS::Iterate)) {
0084     // we dont interpolate, we dont iterate, so we can stop right here.
0085     // cout << "No interpolation, no iteration" << endl;
0086     return *(pairs.begin()->second);
0087   };
0088
0089   // we dont iterate, or we dont have anything to iterate (anymore?)
0090   // so we stop here
0091
0092   // cout << "nobs= " << nobs << "  nq= " << nq << endl;
0093   if (!(theType & SMS::Iterate) || nq <= 2)
0094     return average(pairs, nq);
0095
0096   // we iterate (recursively)
0097
0098   std::vector<GlobalPoint> data1;
0099   std::vector<MyPair>::iterator j;
0100
0101   for (j = pairs.begin(); j - pairs.begin() < nq; ++j)
0102     data1.push_back(*(j->second));
0103
0104   return this->location(data1);
0105 }
0106
0107 GlobalPoint SMS::location(const std::vector<GlPtWt>& wdata) const {
0108   if (!(theType & Weighted)) {
0109     std::vector<GlobalPoint> points;
0110     for (std::vector<GlPtWt>::const_iterator i = wdata.begin(); i != wdata.end(); ++i) {
0111       points.push_back(i->first);
0112     };
0113     if (debug()) {
0114       std::cout << "[SMS] Unweighted SMS was asked for; ignoring the weights." << std::endl;
0115     };
0116     return location(points);
0117   };
0118   // int nobs=wdata.size();
0119   // Sum of weights
0120   float Sumw = 0;
0121   std::vector<GlPtWt>::const_iterator i, j;
0122   for (i = wdata.begin(); i != wdata.end(); ++i)
0123     Sumw += i->second;
0124
0125   // Compute pairwise distances
0126   std::vector<MyPairWt> pairs;
0127   for (i = wdata.begin(); i != wdata.end(); ++i) {
0128     std::vector<FloatPair> D;
0129     // Compute squared distances to all points
0130     for (j = wdata.begin(); j != wdata.end(); ++j)
0131       D.push_back(FloatPair((j->first - i->first).mag2(), j->second));
0132     // Find weighted q-quantile in the distance vector
0133     sort(D.begin(), D.end());
0134     float sumw = 0;
0135     std::vector<FloatPair>::const_iterator where;
0136     for (where = D.begin(); where != D.end(); ++where) {
0137       sumw += where->second;
0138       // cout << sumw << endl;
0139       if (sumw > Sumw * theRatio)
0140         break;
0141     }
0142     MyPairWt tmp(where->first, &(*i));
0143     pairs.push_back(tmp);
0144     // cout << where->first << endl;
0145   };
0146
0147   // Sort pairs by first element
0148   sort(pairs.begin(), pairs.end(), Sorter());
0149
0150   // Find weighted q-quantile in the list of pairs
0151   float sumw = 0;
0152   int nq = 0;
0153   std::vector<MyPairWt>::const_iterator k;
0154   for (k = pairs.begin(); k != pairs.end(); ++k) {
0155     sumw += k->second->second;
0156     ++nq;
0157     if (sumw > Sumw * theRatio)
0158       break;
0159   }
0160
0161   // cout << "nobs= " << nobs << "  nq= " << nq << endl;
0162
0163   if (!(theType & SMS::Interpolate) && !(theType & SMS::Iterate)) {
0164     // we dont interpolate, we dont iterate, so we can stop right here.
0165     // cout << "No interpolation, no iteration" << endl;
0166     return pairs.begin()->second->first;
0167   };
0168
0169   // we dont iterate, or we dont have anything to iterate (anymore?)
0170   // so we stop here
0171
0172   // cout << "nobs= " << nobs << "  nq= " << nq << endl;
0173   if (!(theType & SMS::Iterate) || nq <= 2)
0174     return average(pairs, nq);
0175
0176   // we iterate (recursively)
0177
0178   std::vector<GlPtWt> wdata1;
0179
0180   for (k = pairs.begin(); k - pairs.begin() < nq; ++k)
0181     wdata1.push_back(*(k->second));
0182
0183   return this->location(wdata1);
0184 }