MtvClusterizer1D.h

Go to the documentation of this file.

#ifndef _MtvClusterizer1D_H_
#define _MtvClusterizer1D_H_

#include "CommonTools/Clustering1D/interface/Clusterizer1D.h"
#include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
#include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
#include "CommonTools/Clustering1D/interface/Clustering1DException.h"

#include <vector>
#include <cmath>
#include <algorithm>

template <class T>
class MtvClusterizer1D : public Clusterizer1D<T>
{
public:
    MtvClusterizer1D ( const WeightEstimator<T> & est = TrivialWeightEstimator<T>() );
    MtvClusterizer1D ( const MtvClusterizer1D & );
    ~MtvClusterizer1D();

    std::pair < std::vector < Cluster1D<T> >, std::vector < const T * > >
    operator() ( const std::vector< Cluster1D<T> > & ) const;

    virtual MtvClusterizer1D * clone() const;

private:
    WeightEstimator<T> * theEstimator;
    float theErrorStretchFactor;
};

/*
 *                              --- implementation ---
 *
 */

template <class T>
MtvClusterizer1D<T>::MtvClusterizer1D(
    const MtvClusterizer1D<T> & o ) : theEstimator ( o.theEstimator->clone() )
{}


template <class T>
MtvClusterizer1D<T>::MtvClusterizer1D(
    const WeightEstimator<T> & est ) : theEstimator ( est.clone() )
{}


template <class T>
MtvClusterizer1D<T>::~MtvClusterizer1D()
{
    delete theEstimator;
}

template <class T>
MtvClusterizer1D<T> * MtvClusterizer1D<T>::clone() const
{
    return new MtvClusterizer1D<T> ( * this );
}

template <class T>
std::pair < std::vector < Cluster1D<T> >, std::vector < const T * > >
MtvClusterizer1D<T>::operator() ( const std::vector < Cluster1D<T> > & ov ) const
{
    typedef Cluster1D<T> Cluster1D;
    using namespace Clusterizer1DCommons;
    std::vector < const T * > unusedtracks;
    switch ( ov.size() )
    {
    case 0:
        throw Clustering1DException("[MtvClusterizer1D] no values given" );
    case 1:
        std::pair < std::vector < Cluster1D >, std::vector < const T * > > ret ( ov, unusedtracks );
        return ret;
    };
    std::vector < Cluster1D > v = ov;
    sort ( v.begin(), v.end(), ComparePairs<T>() );
    std::vector < Cluster1D > sols;
    std::vector < const T * > trks;

    typename std::vector< Cluster1D >::iterator cur = v.begin();
    typename std::vector< Cluster1D >::iterator   end = (v.end() - 1 );
    double cur_min = cur->weight() + ( cur+1 )->weight();

    for ( typename std::vector< Cluster1D >::iterator i=v.begin();
            i!=end ; ++i )
    {
        double cur_val = i->weight() + ( i+1 )->weight();
        if ( cur_val > cur_min )
        {
            cur_min = cur_val;
            cur = i;
        };
    };

    double weight = ( cur->weight() + (cur+1)->weight() );
    double est = ( cur->weight() * cur->position().value() +
                   (cur+1)->weight() * (cur+1)->position().value()) / weight;
    double sigma = sqrt ( square ( cur->position().value() - est ) +
                          square ( (cur+1)->position().value() - est ) );
    double err=0.;
    int inliers=0;

    for ( typename std::vector< Cluster1D >::iterator i=v.begin();
            i!=v.end() ; ++i )
    {
        if ( fabs ( i->position().value() - est ) < 3 * sigma )
        {
            // all within 3 sigma are 'in'
            add
                ( i->tracks(), trks );
            err+= square ( i->position().value() - est );
            inliers++;
        }
        else
        {
            add
                ( i->tracks(), unusedtracks );
        };
    };
    err /= ( inliers - 1 ); // the algo definitely produces 2 or more inliers
    err = sqrt ( err );

    sols.push_back ( Cluster1D ( Measurement1D ( est,err ), trks, weight ) );
    std::pair < std::vector < Cluster1D >, std::vector < const T * > > ret ( sols, unusedtracks );
    return ret;
}

#endif