PixelVertexFinding/interface/DivisiveClusterizer1D.h

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