Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:26:51

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