Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-06-12 03:09:04

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   int ncount = 0;
0106   std::vector<Cluster1D<T> > partOfPTracks;
0107   partOfPTracks.push_back(input.front());
0108   for (typename std::vector<Cluster1D<T> >::const_iterator ic = (input.begin()) + 1; ic != input.end(); ic++) {
0109     ncount++;
0110     if (fabs((*ic).position().value() - (*(ic - 1)).position().value()) < (double)theZSeparation) {
0111       partOfPTracks.push_back((*ic));
0112     } else {
0113       if (partOfPTracks.size() >= theNTkMin) {
0114         std::vector<Cluster1D<T> > clusters = makeCluster1Ds(partOfPTracks, discardedTracks);
0115         for (typename std::vector<Cluster1D<T> >::const_iterator iclus = clusters.begin(); iclus != clusters.end();
0116              iclus++) {
0117           finalCluster1Ds.push_back(*iclus);
0118         }
0119         insertTracks(discardedTracks, totDiscardedTracks);
0120       } else {
0121         insertTracks(partOfPTracks, totDiscardedTracks);
0122       }
0123       partOfPTracks.clear();
0124       partOfPTracks.push_back((*ic));
0125     }
0126   }
0127   if (partOfPTracks.size() >= theNTkMin) {
0128     std::vector<Cluster1D<T> > clusters = makeCluster1Ds(partOfPTracks, discardedTracks);
0129     for (typename std::vector<Cluster1D<T> >::const_iterator iclus = clusters.begin(); iclus != clusters.end();
0130          iclus++) {
0131       finalCluster1Ds.push_back(*iclus);
0132     }
0133     insertTracks(discardedTracks, totDiscardedTracks);
0134   } else {
0135     insertTracks(partOfPTracks, totDiscardedTracks);
0136   }
0137 
0138   sort(finalCluster1Ds.begin(), finalCluster1Ds.end(), ComparePairs<T>());
0139   // reverse(theCluster1Ds.begin(), theCluster1Ds.end());
0140 
0141   return;
0142 }
0143 
0144 template <class T>
0145 std::vector<Cluster1D<T> > DivisiveClusterizer1D<T>::makeCluster1Ds(std::vector<Cluster1D<T> >& clusters,
0146                                                                     std::vector<Cluster1D<T> >& discardedTracks) const {
0147   std::vector<Cluster1D<T> > finalCluster1Ds;
0148   discardedTracks.clear();
0149   std::vector<Cluster1D<T> > pvClu0 = clusters;
0150   std::vector<Cluster1D<T> > pvCluNew = pvClu0;
0151   bool stop = false;
0152   while (!stop) {
0153     int nDiscardedAtIteration = 100;
0154     while (nDiscardedAtIteration != 0) {
0155       pvCluNew = theCleaner->clusters(pvClu0);
0156       std::vector<Cluster1D<T> > tracksAtIteration = theCleaner->discardedCluster1Ds();
0157       nDiscardedAtIteration = tracksAtIteration.size();
0158       if (nDiscardedAtIteration != 0) {
0159         insertTracks(tracksAtIteration, discardedTracks);
0160         pvClu0 = pvCluNew;
0161       }
0162     }  // while nDiscardedAtIteration
0163     unsigned ntkclus = pvCluNew.size();
0164     unsigned ndiscard = discardedTracks.size();
0165 
0166     if (ntkclus >= theNTkMin) {
0167       //save the cluster
0168       finalCluster1Ds.push_back(mergeCluster1Ds(pvCluNew));
0169       if (ndiscard >= theNTkMin) {  //make a new cluster and reset
0170         pvClu0 = discardedTracks;
0171         discardedTracks.clear();
0172       } else {  //out of loop
0173         stop = true;
0174       }
0175     } else {
0176       insertTracks(pvCluNew, discardedTracks);
0177       stop = true;
0178     }
0179   }  // while stop
0180   return finalCluster1Ds;
0181 }
0182 
0183 template <class T>
0184 void DivisiveClusterizer1D<T>::insertTracks(std::vector<Cluster1D<T> >& clusou,
0185                                             std::vector<Cluster1D<T> >& cludest) const {
0186   if (clusou.empty())
0187     return;
0188   for (typename std::vector<Cluster1D<T> >::const_iterator iclu = clusou.begin(); iclu != clusou.end(); iclu++) {
0189     cludest.push_back(*iclu);
0190   }
0191   /*
0192     for ( typename std::vector< Cluster1D<T> >::const_iterator iclu = clu.begin(); 
0193     iclu != clu.end(); iclu++){
0194       if (total) {
0195         theTotalDiscardedTracks.push_back(*iclu);
0196       }else { 
0197         theDiscardedTracks.push_back(*iclu);
0198       }
0199     }
0200     */
0201   return;
0202 }
0203 
0204 template <class T>
0205 std::vector<const T*> DivisiveClusterizer1D<T>::takeTracks(const std::vector<Cluster1D<T> >& clu) const {
0206   std::vector<const T*> tracks;
0207   for (typename std::vector<Cluster1D<T> >::const_iterator iclu = clu.begin(); iclu != clu.end(); iclu++) {
0208     std::vector<const T*> clutks = iclu->tracks();
0209     for (typename std::vector<const T*>::const_iterator i = clutks.begin(); i != clutks.end(); ++i) {
0210       tracks.push_back(*i);
0211     }
0212   }
0213   return tracks;
0214 }
0215 
0216 template <class T>
0217 Cluster1D<T> DivisiveClusterizer1D<T>::mergeCluster1Ds(std::vector<Cluster1D<T> >& clusters) const {
0218   Cluster1D<T> result = clusters.front();
0219   for (typename std::vector<Cluster1D<T> >::iterator iclu = (clusters.begin()) + 1; iclu != clusters.end(); iclu++) {
0220     Cluster1D<T> old = result;
0221     result = (*theMerger)(old, *iclu);
0222   }
0223   return result;
0224 }
0225 
0226 #endif