FuzzyKMeanAlg.h

Go to the documentation of this file.

///
/// \file    FuzzyKMeanAlg.h
/// \brief   Utility class for making prongs from vertex
/// \version $Id: FuzzyKMeanAlg.h,v 1.4 2012-10-01 00:29:02 edniner Exp $
/// \author  eniner@indiana.edu
///
#ifndef FUZZ_FUZZYKMEANALG_H
#define FUZZ_FUZZYKMEANALG_H

#include "fhiclcpp/types/Atom.h"

#include "RecoBase/HitList.h"

#include <string>

namespace rb
{
  class Cluster;
  class HitList;
}


namespace fuzz
{
  struct FuzzyKMeanParams
  {
    template<class T> using Atom = fhicl::Atom<T>;
    using Comment = fhicl::Comment;
    using Name = fhicl::Name;

    Atom<bool> WeightByQ{Name("WeightByQ"), Comment("Weight all hits equally (false), or weight by charge (true)")};
    Atom<double> Fuzz{Name("Fuzz"), Comment("Degree of fuzzykness in clusters")};
    Atom<double> Tol{Name("Tol"), Comment("Tolerance for convergence of cluster centers on a solution")};
    Atom<double> Beta{Name("Beta"), Comment("Cluster spread (2 sigma)")};
    Atom<double> Memb{Name("Memb"), Comment("Min membership to a clust")};
    Atom<unsigned int> MaxItr{Name("Maxitr"), Comment("Maximum number of clustering iterations")};

    Atom<double> Par1{Name("Par1"), Comment("Angular uncertainty has the form F(x)=Par1/x + Par2 + Par3*x."
                                            "The parameters were modeled based off of simulations of muons and electrons.")};
    Atom<double> Par2{Name("Par2"), Comment("Angular uncertainty has the form F(x)=Par1/x + Par2 + Par3*x."
                                            "The parameters were modeled based off of simulations of muons and electrons.")};
    Atom<double> Par3{Name("Par3"), Comment("Angular uncertainty has the form F(x)=Par1/x + Par2 + Par3*x."
                                            "The parameters were modeled based off of simulations of muons and electrons.")};

    Atom<int> Bins{Name("Bins"), Comment("Number of angular bins to use to make density map for seeding")};
    Atom<double> MergeAng{Name("MergeAng"), Comment("Maximum angular difference between cluster centers to consider a merge (radians)")};
    Atom<double> MergeMemb{Name("MergeMemb"), Comment("Minimum degree of sharing between two clusters in order to merge")};

    Atom<double> Cutoff{Name("CutOff"), Comment("Distance (in cm) beyond which the angular uncertainty should remain constant")};
    Atom<double> Duplicate{Name("Duplicate"), Comment("Threshold for two cluster centers being duplicates (in radians)")};
    Atom<int> PrintLevel{Name("PrintLevel"), Comment("Level of debugging to control print statements (0=no debug, 1= major info, 2= minor info)")};
  };


  class FuzzyKMeanAlg
  {
  public:
    explicit FuzzyKMeanAlg(const FuzzyKMeanParams& params);
    ~FuzzyKMeanAlg();

    /// \brief Function to take in a cluster and build the hitlist needed to do the clustering
    void LoadCluster(rb::Cluster slice);

    /// \brief A function to control all steps in clustering in 2D on a group of hits
    /// \param minClust The minimum number of cluster candidates to make
    /// \param maxClust The maximum number of cluster candidates to make
    /// \param xy The x or y coordinate (depending on view) of the vertex to cluster around
    /// \param z The z coordinate of the vertex to cluster around
    /// \param view The detector view (either kX or kY) of the hits to use in clustering. This must be a 2d view
    void   DoClustering(unsigned int minClust, unsigned int maxClust, 
                        double xy, double z, geo::View_t view);

    /// \brief A function to calculate membership weights for a specified number of clusters
    /// \param n The number of cluster centers to use 
    void   MakeClusters(unsigned int n);

    /// \brief A function to calculate the seeds used for the first round of clustering
    /// \param n Number of seeds to calculate
    void   SeedA(unsigned int n);

    /// \brief A function to collect the hits that have not made it to a cluster yet and assign new cluster seeds based on that
    bool   FindUnclustered();

    /// \brief Calculate and return the distance from a hit to a cluster center
    /// \param mu Mean position of cluster center (in radians)
    /// \param x Position of a hit (in radians)
    /// \param dx The uncertainty in hit position (in radians)
    double Distance(double mu, double x, double dx);

    /// \brief Compute all distances for j hits to i cluster centers
    void   ComputeDij();

    /// \brief Compute the membership values for j hits to i cluster centers
    void   ComputeUij();

    /// \brief Update the mean cluster center positions
    bool   ComputeA();

    /// \brief Calculate and return the score evaluating the clustering
    double Jpfuz();

    /// \brief Check if each of the cluster centers (after the algorith has converged) is unique, remove duplicates to reseed
    bool   CheckUnique();

    /// \brief Function to turn the 2D hit coordinates into angles (radians) with respect to the vertex for clustering
    /// \param vXY The x or y (depending on view) vertex coordinate
    /// \param vZ The z vertex coordinate
    /// \param view The view, either kX or kY, of the hits to use.  This needs to match the view supplied in vXY
    void   MakeAngles(double vXY, double vZ, geo::View_t view);

    /// \brief Function to merge cluster centers together that have significant overlap
    void   MergeClusters();

  public:
    /// List of hits in slice
    rb::HitList                            fHL;    

    /// Distance from vertex
    std::vector<double>                    fDist;  

    /// Hit weight, currently PE
    std::vector<double>                    fW;   

    /// Angle to vertex  
    std::vector<double>                    fX;    

    /// Uncertainty in angle 
    std::vector<double>                    fdX;    

    /// Dist clust i to sub. j
    std::vector< std::vector<double> >     fD;     

    /// Assign cluster i to sub. j
    std::vector< std::vector<double> >     fU;     

    /// Cluster centers
    std::vector<double>                    fA;     
  
    /// List of used seeds
    std::vector<double>                    fSeed;  

    /// Unused cell pos., uncert
    std::vector<std::pair<double,double> > fUnused;

    /// List of max densities
    std::vector<std::pair<double,double> > fLMax;  

  private:
    FuzzyKMeanParams fParams;
  };
}


#endif
////////////////////////////////////////////////////////////////////////