Clustering1D/interface/MtvClusterizer1D.h

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 /** Hsm clusterizer in one dimension, originally designed for ApexPoint Finding
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  *                              --- implementation ---
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       // all within 3 sigma are 'in'
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);  // the algo definitely produces 2 or more inliers
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