PixelVertexFinding/interface/Cluster1DCleaner.h

0001 #ifndef RecoVertex_PixelVertexFinding_Cluster1DCleaner_h
0002 #define RecoVertex_PixelVertexFinding_Cluster1DCleaner_h
0003
0004 #include "CommonTools/Clustering1D/interface/Cluster1D.h"
0005
0006 #include <cmath>
0007
0008 /*
0009  * given a vector<Cluster1D<T> >, erase Cluster1D further away than
0010  * ZOffeSet from the average position, then
0011  * recompute the vertex position. The ZOffeSet is taken as
0012  * an aurgument
0013 */
0014 namespace pixeltemp {
0015   template <class T>
0016   class Cluster1DCleaner {
0017   public:
0018     Cluster1DCleaner(const float zoffset, bool useErr) : theZOffSet(zoffset), theUseError(useErr) {
0019       theCleanedCluster1Ds.clear();
0020       theDiscardedCluster1Ds.clear();
0021     }
0022     // return the compatible clusters
0023     std::vector<Cluster1D<T> > clusters(const std::vector<Cluster1D<T> >&);
0024     /*
0025        return the vector of discarded Cluster1Ds
0026        it should be called after Cluster1DCleaner::clusters
0027        otherwise return an empty vector
0028     */
0029     std::vector<Cluster1D<T> > discardedCluster1Ds() const { return theDiscardedCluster1Ds; }
0030
0031   private:
0032     void cleanCluster1Ds(const std::vector<Cluster1D<T> >&);
0033     float average(const std::vector<Cluster1D<T> >&);
0034     std::vector<Cluster1D<T> > theCleanedCluster1Ds;
0035     std::vector<Cluster1D<T> > theDiscardedCluster1Ds;
0036     float theZOffSet;
0037     bool theUseError;
0038   };
0039
0040   /*
0041  *                                implementation
0042  */
0043
0044   template <class T>
0045   std::vector<Cluster1D<T> > Cluster1DCleaner<T>::clusters(const std::vector<Cluster1D<T> >& clust) {
0046     cleanCluster1Ds(clust);
0047     return theCleanedCluster1Ds;
0048   }
0049
0050   template <class T>
0051   void Cluster1DCleaner<T>::cleanCluster1Ds(const std::vector<Cluster1D<T> >& clust) {
0052     theCleanedCluster1Ds.clear();
0053     theDiscardedCluster1Ds.clear();
0054     if (clust.empty())
0055       return;
0056     float oldPos = average(clust);
0057     for (typename std::vector<Cluster1D<T> >::const_iterator ic = clust.begin(); ic != clust.end(); ic++) {
0058       float discr = theUseError ? fabs(((*ic).position().value() - oldPos) / (*ic).position().error())
0059                                 : fabs(((*ic).position().value() - oldPos));
0060       if (discr < theZOffSet) {
0061         theCleanedCluster1Ds.push_back(*ic);
0062       } else {
0063         theDiscardedCluster1Ds.push_back(*ic);
0064       }
0065     }
0066     return;
0067   }
0068
0069   //I could use the Cluster1DMerger...
0070   template <class T>
0071   float Cluster1DCleaner<T>::average(const std::vector<Cluster1D<T> >& clust) {
0072     //    float ave = clust.front().position().value();
0073     //    float err = clust.front().position().error();
0074     //    for( typename std::vector < Cluster1D<T> >::const_iterator ic=(clust.begin())+1;
0075     //         ic != clust.end(); ic++)
0076     //    {
0077     //        float oldave = ave;
0078     //        float olderr = err;
0079     //        ave = ( oldave/olderr/olderr +
0080     //                ic->position().value()/ic->position().error()/ic->position().error()) /
0081     //              (1./olderr/olderr + 1./ic->position().error()/ic->position().error());
0082     //        err = sqrt(olderr*olderr + ic->position().error()*ic->position().error());
0083     //    }
0084     float sumUp = 0;
0085     float sumDown = 0;
0086     for (typename std::vector<Cluster1D<T> >::const_iterator ic = (clust.begin()) + 1; ic != clust.end(); ic++) {
0087       float err2 = ic->position().error();
0088       err2 *= err2;
0089       if (err2 != 0) {
0090         sumUp += ic->position().value() / err2;  // error-weighted average of Z at IP
0091         sumDown += 1 / err2;
0092       }
0093     }
0094     return (sumDown > 0) ? sumUp / sumDown : 0;
0095   }
0096 }  // namespace pixeltemp
0097 #endif