Clustering1D/interface/DivisiveClusterizer1D.h

0001 #ifndef _DivisiveClusterizer1D_H_
0002 #define _DivisiveClusterizer1D_H_
0003
0004 #include "CommonTools/Clustering1D/interface/Clusterizer1D.h"
0005 #include "CommonTools/Clustering1D/interface/Cluster1DMerger.h"
0006 #include "CommonTools/Clustering1D/interface/Cluster1DCleaner.h"
0007 #include "CommonTools/Clustering1D/interface/TrivialWeightEstimator.h"
0008 #include "CommonTools/Clustering1D/interface/Clusterizer1DCommons.h"
0009
0010 /**
0011  * Find the modes with a simple divisive method.
0012  */
0013 template <class T>
0014 class DivisiveClusterizer1D : public Clusterizer1D<T> {
0015 public:
0016   /**
0017      *  \param zoffset maximum distance between track position and position of its cluster
0018      *         (depending on useError its either weighted or physical distance)
0019      *  \param ntkmin Minimum number of tracks required to form a cluster.
0020      *  \param useError physical distances or weighted distances.
0021      *  \param zsep Maximum distance between two adjacent tracks that belong
0022      *         to the same initial cluster.
0023      *  \param wei Compute the cluster "center" with an unweighted or a weighted
0024      *         average of the tracks. Weighted means weighted with the error
0025      *         of the data point.
0026      */
0027   DivisiveClusterizer1D(float zoffset = 5., int ntkmin = 5, bool useError = true, float zsep = 0.05, bool wei = true);
0028
0029   ~DivisiveClusterizer1D();
0030
0031   std::pair<std::vector<Cluster1D<T> >, std::vector<const T*> > operator()(const std::vector<Cluster1D<T> >&) const;
0032   virtual DivisiveClusterizer1D* clone() const;
0033
0034 private:
0035   //methods
0036   void findCandidates(const std::vector<Cluster1D<T> >&,
0037                       std::vector<Cluster1D<T> >&,
0038                       std::vector<Cluster1D<T> >&) const;
0039   std::vector<Cluster1D<T> > makeCluster1Ds(std::vector<Cluster1D<T> >&, std::vector<Cluster1D<T> >&) const;
0040   void insertTracks(std::vector<Cluster1D<T> >&, std::vector<Cluster1D<T> >&) const;
0041   std::vector<const T*> takeTracks(const std::vector<Cluster1D<T> >&) const;
0042   Cluster1D<T> mergeCluster1Ds(std::vector<Cluster1D<T> >&) const;
0043   //data members
0044   // std::vector<Cluster1D<T> > theCluster1Ds;
0045   // std::vector<Cluster1D<T> > theTotalDiscardedTracks;
0046   //  std::vector<Cluster1D<T> > theDiscardedTracks;
0047   Cluster1DMerger<T>* theMerger;
0048   Cluster1DCleaner<T>* theCleaner;
0049   float theZOffSet;
0050   float theZSeparation;
0051   unsigned theNTkMin;
0052   bool theWei;
0053   bool theUseError;
0054 };
0055
0056 /*
0057  *        implementation
0058  *
0059  */
0060
0061 template <class T>
0062 DivisiveClusterizer1D<T>::DivisiveClusterizer1D(float zoffset, int ntkmin, bool useError, float zsep, bool wei)
0063     : theZOffSet(zoffset), theZSeparation(zsep), theNTkMin(ntkmin), theWei(wei), theUseError(useError) {
0064   //  theDiscardedTracks.clear();
0065   // theTotalDiscardedTracks.clear();
0066   //  theCluster1Ds.clear();
0067   TrivialWeightEstimator<T> weightEstimator;
0068   theMerger = new Cluster1DMerger<T>(weightEstimator);
0069   theCleaner = new Cluster1DCleaner<T>(theZOffSet, theUseError);
0070 }
0071
0072 template <class T>
0073 DivisiveClusterizer1D<T>::~DivisiveClusterizer1D() {
0074   delete theMerger;
0075   delete theCleaner;
0076 }
0077
0078 template <class T>
0079 std::pair<std::vector<Cluster1D<T> >, std::vector<const T*> > DivisiveClusterizer1D<T>::operator()(
0080     const std::vector<Cluster1D<T> >& input) const {
0081   std::vector<Cluster1D<T> > discardedCluster1Ds;
0082   std::vector<Cluster1D<T> > output;
0083   findCandidates(input, output, discardedCluster1Ds);
0084   return std::pair<std::vector<Cluster1D<T> >, std::vector<const T*> >(output, takeTracks(discardedCluster1Ds));
0085 }
0086
0087 template <class T>
0088 DivisiveClusterizer1D<T>* DivisiveClusterizer1D<T>::clone() const {
0089   return new DivisiveClusterizer1D<T>(*this);
0090 }
0091
0092 template <class T>
0093 void DivisiveClusterizer1D<T>::findCandidates(const std::vector<Cluster1D<T> >& inputo,
0094                                               std::vector<Cluster1D<T> >& finalCluster1Ds,
0095                                               std::vector<Cluster1D<T> >& totDiscardedTracks) const {
0096   using namespace Clusterizer1DCommons;
0097
0098   std::vector<Cluster1D<T> > input = inputo;
0099   std::vector<Cluster1D<T> > discardedTracks;
0100   if (input.size() < theNTkMin) {
0101     insertTracks(input, totDiscardedTracks);
0102     return;
0103   }
0104   sort(input.begin(), input.end(), ComparePairs<T>());
0105   std::vector<Cluster1D<T> > partOfPTracks;
0106   partOfPTracks.push_back(input.front());
0107   for (typename std::vector<Cluster1D<T> >::const_iterator ic = (input.begin()) + 1; ic != input.end(); ic++) {
0108     if (fabs((*ic).position().value() - (*(ic - 1)).position().value()) < (double)theZSeparation) {
0109       partOfPTracks.push_back((*ic));
0110     } else {
0111       if (partOfPTracks.size() >= theNTkMin) {
0112         std::vector<Cluster1D<T> > clusters = makeCluster1Ds(partOfPTracks, discardedTracks);
0113         for (typename std::vector<Cluster1D<T> >::const_iterator iclus = clusters.begin(); iclus != clusters.end();
0114              iclus++) {
0115           finalCluster1Ds.push_back(*iclus);
0116         }
0117         insertTracks(discardedTracks, totDiscardedTracks);
0118       } else {
0119         insertTracks(partOfPTracks, totDiscardedTracks);
0120       }
0121       partOfPTracks.clear();
0122       partOfPTracks.push_back((*ic));
0123     }
0124   }
0125   if (partOfPTracks.size() >= theNTkMin) {
0126     std::vector<Cluster1D<T> > clusters = makeCluster1Ds(partOfPTracks, discardedTracks);
0127     for (typename std::vector<Cluster1D<T> >::const_iterator iclus = clusters.begin(); iclus != clusters.end();
0128          iclus++) {
0129       finalCluster1Ds.push_back(*iclus);
0130     }
0131     insertTracks(discardedTracks, totDiscardedTracks);
0132   } else {
0133     insertTracks(partOfPTracks, totDiscardedTracks);
0134   }
0135
0136   sort(finalCluster1Ds.begin(), finalCluster1Ds.end(), ComparePairs<T>());
0137   // reverse(theCluster1Ds.begin(), theCluster1Ds.end());
0138
0139   return;
0140 }
0141
0142 template <class T>
0143 std::vector<Cluster1D<T> > DivisiveClusterizer1D<T>::makeCluster1Ds(std::vector<Cluster1D<T> >& clusters,
0144                                                                     std::vector<Cluster1D<T> >& discardedTracks) const {
0145   std::vector<Cluster1D<T> > finalCluster1Ds;
0146   discardedTracks.clear();
0147   std::vector<Cluster1D<T> > pvClu0 = clusters;
0148   std::vector<Cluster1D<T> > pvCluNew = pvClu0;
0149   bool stop = false;
0150   while (!stop) {
0151     int nDiscardedAtIteration = 100;
0152     while (nDiscardedAtIteration != 0) {
0153       pvCluNew = theCleaner->clusters(pvClu0);
0154       std::vector<Cluster1D<T> > tracksAtIteration = theCleaner->discardedCluster1Ds();
0155       nDiscardedAtIteration = tracksAtIteration.size();
0156       if (nDiscardedAtIteration != 0) {
0157         insertTracks(tracksAtIteration, discardedTracks);
0158         pvClu0 = pvCluNew;
0159       }
0160     }  // while nDiscardedAtIteration
0161     unsigned ntkclus = pvCluNew.size();
0162     unsigned ndiscard = discardedTracks.size();
0163
0164     if (ntkclus >= theNTkMin) {
0165       //save the cluster
0166       finalCluster1Ds.push_back(mergeCluster1Ds(pvCluNew));
0167       if (ndiscard >= theNTkMin) {  //make a new cluster and reset
0168         pvClu0 = discardedTracks;
0169         discardedTracks.clear();
0170       } else {  //out of loop
0171         stop = true;
0172       }
0173     } else {
0174       insertTracks(pvCluNew, discardedTracks);
0175       stop = true;
0176     }
0177   }  // while stop
0178   return finalCluster1Ds;
0179 }
0180
0181 template <class T>
0182 void DivisiveClusterizer1D<T>::insertTracks(std::vector<Cluster1D<T> >& clusou,
0183                                             std::vector<Cluster1D<T> >& cludest) const {
0184   if (clusou.empty())
0185     return;
0186   for (typename std::vector<Cluster1D<T> >::const_iterator iclu = clusou.begin(); iclu != clusou.end(); iclu++) {
0187     cludest.push_back(*iclu);
0188   }
0189   /*
0190     for ( typename std::vector< Cluster1D<T> >::const_iterator iclu = clu.begin();
0191     iclu != clu.end(); iclu++){
0192       if (total) {
0193         theTotalDiscardedTracks.push_back(*iclu);
0194       }else {
0195         theDiscardedTracks.push_back(*iclu);
0196       }
0197     }
0198     */
0199   return;
0200 }
0201
0202 template <class T>
0203 std::vector<const T*> DivisiveClusterizer1D<T>::takeTracks(const std::vector<Cluster1D<T> >& clu) const {
0204   std::vector<const T*> tracks;
0205   for (typename std::vector<Cluster1D<T> >::const_iterator iclu = clu.begin(); iclu != clu.end(); iclu++) {
0206     std::vector<const T*> clutks = iclu->tracks();
0207     for (typename std::vector<const T*>::const_iterator i = clutks.begin(); i != clutks.end(); ++i) {
0208       tracks.push_back(*i);
0209     }
0210   }
0211   return tracks;
0212 }
0213
0214 template <class T>
0215 Cluster1D<T> DivisiveClusterizer1D<T>::mergeCluster1Ds(std::vector<Cluster1D<T> >& clusters) const {
0216   Cluster1D<T> result = clusters.front();
0217   for (typename std::vector<Cluster1D<T> >::iterator iclu = (clusters.begin()) + 1; iclu != clusters.end(); iclu++) {
0218     Cluster1D<T> old = result;
0219     result = (*theMerger)(old, *iclu);
0220   }
0221   return result;
0222 }
0223
0224 #endif