Back to home page

Project CMSSW displayed by LXR

 
 

    


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 /** 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