File indexing completed on 2023-03-17 11:23:35
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 }
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
0065
0066
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
0072 for (std::vector<GlobalPoint>::const_iterator j = data.begin(); j != data.end(); ++j) {
0073 D.push_back((*j - *i).mag2());
0074 }
0075
0076 sort(D.begin(), D.end());
0077 MyPair tmp(D[nq - 1], &(*i));
0078 pairs.push_back(tmp);
0079 };
0080
0081
0082 sort(pairs.begin(), pairs.end(), Sorter());
0083 if (!(theType & SMS::Interpolate) && !(theType & SMS::Iterate)) {
0084
0085
0086 return *(pairs.begin()->second);
0087 };
0088
0089
0090
0091
0092
0093 if (!(theType & SMS::Iterate) || nq <= 2)
0094 return average(pairs, nq);
0095
0096
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
0119
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
0126 std::vector<MyPairWt> pairs;
0127 for (i = wdata.begin(); i != wdata.end(); ++i) {
0128 std::vector<FloatPair> D;
0129
0130 for (j = wdata.begin(); j != wdata.end(); ++j)
0131 D.push_back(FloatPair((j->first - i->first).mag2(), j->second));
0132
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
0139 if (sumw > Sumw * theRatio)
0140 break;
0141 }
0142 MyPairWt tmp(where->first, &(*i));
0143 pairs.push_back(tmp);
0144
0145 };
0146
0147
0148 sort(pairs.begin(), pairs.end(), Sorter());
0149
0150
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
0162
0163 if (!(theType & SMS::Interpolate) && !(theType & SMS::Iterate)) {
0164
0165
0166 return pairs.begin()->second->first;
0167 };
0168
0169
0170
0171
0172
0173 if (!(theType & SMS::Iterate) || nq <= 2)
0174 return average(pairs, nq);
0175
0176
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 }