File indexing completed on 2024-04-06 12:01:04
0001 #ifndef _MtvClusterizer1D_H_
0002 #define _MtvClusterizer1D_H_
0003
0004 #include "CommonTools/Clustering1D/interface/Clusterizer1D.h"
0005 #include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
0006 #include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
0007 #include "CommonTools/Clustering1D/interface/Clustering1DException.h"
0008
0009 #include <vector>
0010 #include <cmath>
0011 #include <algorithm>
0012
0013
0014
0015
0016 template <class T>
0017 class MtvClusterizer1D : public Clusterizer1D<T> {
0018 public:
0019 MtvClusterizer1D(const WeightEstimator<T>& est = TrivialWeightEstimator<T>());
0020 MtvClusterizer1D(const MtvClusterizer1D&);
0021 ~MtvClusterizer1D();
0022
0023 std::pair<std::vector<Cluster1D<T> >, std::vector<const T*> > operator()(const std::vector<Cluster1D<T> >&) const;
0024
0025 virtual MtvClusterizer1D* clone() const;
0026
0027 private:
0028 WeightEstimator<T>* theEstimator;
0029 float theErrorStretchFactor;
0030 };
0031
0032
0033
0034
0035
0036
0037 template <class T>
0038 MtvClusterizer1D<T>::MtvClusterizer1D(const MtvClusterizer1D<T>& o) : theEstimator(o.theEstimator->clone()) {}
0039
0040 template <class T>
0041 MtvClusterizer1D<T>::MtvClusterizer1D(const WeightEstimator<T>& est) : theEstimator(est.clone()) {}
0042
0043 template <class T>
0044 MtvClusterizer1D<T>::~MtvClusterizer1D() {
0045 delete theEstimator;
0046 }
0047
0048 template <class T>
0049 MtvClusterizer1D<T>* MtvClusterizer1D<T>::clone() const {
0050 return new MtvClusterizer1D<T>(*this);
0051 }
0052
0053 template <class T>
0054 std::pair<std::vector<Cluster1D<T> >, std::vector<const T*> > MtvClusterizer1D<T>::operator()(
0055 const std::vector<Cluster1D<T> >& ov) const {
0056 typedef Cluster1D<T> Cluster1D;
0057 using namespace Clusterizer1DCommons;
0058 std::vector<const T*> unusedtracks;
0059 switch (ov.size()) {
0060 case 0:
0061 throw Clustering1DException("[MtvClusterizer1D] no values given");
0062 case 1:
0063 std::pair<std::vector<Cluster1D>, std::vector<const T*> > ret(ov, unusedtracks);
0064 return ret;
0065 };
0066 std::vector<Cluster1D> v = ov;
0067 sort(v.begin(), v.end(), ComparePairs<T>());
0068 std::vector<Cluster1D> sols;
0069 std::vector<const T*> trks;
0070
0071 typename std::vector<Cluster1D>::iterator cur = v.begin();
0072 typename std::vector<Cluster1D>::iterator end = (v.end() - 1);
0073 double cur_min = cur->weight() + (cur + 1)->weight();
0074
0075 for (typename std::vector<Cluster1D>::iterator i = v.begin(); i != end; ++i) {
0076 double cur_val = i->weight() + (i + 1)->weight();
0077 if (cur_val > cur_min) {
0078 cur_min = cur_val;
0079 cur = i;
0080 };
0081 };
0082
0083 double weight = (cur->weight() + (cur + 1)->weight());
0084 double est = (cur->weight() * cur->position().value() + (cur + 1)->weight() * (cur + 1)->position().value()) / weight;
0085 double sigma = sqrt(square(cur->position().value() - est) + square((cur + 1)->position().value() - est));
0086 double err = 0.;
0087 int inliers = 0;
0088
0089 for (typename std::vector<Cluster1D>::iterator i = v.begin(); i != v.end(); ++i) {
0090 if (fabs(i->position().value() - est) < 3 * sigma) {
0091
0092 add(i->tracks(), trks);
0093 err += square(i->position().value() - est);
0094 inliers++;
0095 } else {
0096 add(i->tracks(), unusedtracks);
0097 };
0098 };
0099 err /= (inliers - 1);
0100 err = sqrt(err);
0101
0102 sols.push_back(Cluster1D(Measurement1D(est, err), trks, weight));
0103 std::pair<std::vector<Cluster1D>, std::vector<const T*> > ret(sols, unusedtracks);
0104 return ret;
0105 }
0106
0107 #endif