PixelVertexFinding/interface/FindPeakFastPV.h

0001 #ifndef RecoVertex_PixelVertexFinding_FindPeakFastPV_h
0002 #define RecoVertex_PixelVertexFinding_FindPeakFastPV_h
0003 /** \class FindPeakFastPV FindPeakFastPV.h RecoVertex/PixelVertexFinding/FindPeakFastPV.h
0004  * Given *zProjections* and *zWeights* find the peak of width *m_zClusterWidth*.
0005  * Use only values with *zWeights*>*m_weightCut*.
0006  * Look near *oldVertex* within *m_zClusterSearchArea*.
0007  *
0008  *  \author Silvio Donato (SNS)
0009  */
0010
0011 #include <vector>
0012
0013 inline float FindPeakFastPV(const std::vector<float> &zProjections,
0014                             const std::vector<float> &zWeights,
0015                             const float oldVertex,
0016                             const float m_zClusterWidth,
0017                             const float m_zClusterSearchArea,
0018                             const float m_weightCut) {
0019   float centerWMax = oldVertex;
0020   if (m_zClusterWidth > 0 && m_zClusterSearchArea > 0) {
0021     std::vector<float>::const_iterator itCenter = zProjections.begin();
0022     std::vector<float>::const_iterator itLeftSide = zProjections.begin();
0023     std::vector<float>::const_iterator itRightSide = zProjections.begin();
0024     const float zClusterWidth = m_zClusterWidth * 0.5;  //take half zCluster width
0025     const float zLowerBound = oldVertex - zClusterWidth - m_zClusterSearchArea;
0026     const float zUpperBound = oldVertex + zClusterWidth + m_zClusterSearchArea;
0027     float maxW = 0;
0028     for (; itCenter != zProjections.end(); itCenter++) {
0029       //loop only on the zProjections within oldVertex +/- (zClusterWidth + m_zClusterSearchArea)
0030       if ((*itCenter < zLowerBound) || (*itCenter > zUpperBound))
0031         continue;
0032
0033       while (itLeftSide != zProjections.end() && (*itCenter - *itLeftSide) > zClusterWidth)
0034         itLeftSide++;
0035       while (itRightSide != zProjections.end() && (*itRightSide - *itCenter) < zClusterWidth)
0036         itRightSide++;
0037       float nWeighted = 0;
0038       float centerW = 0;
0039
0040       for (std::vector<float>::const_iterator ii = itLeftSide; ii != itRightSide; ii++) {
0041         //loop inside the peak and calculate its weight
0042         if (zWeights[ii - zProjections.begin()] < m_weightCut)
0043           continue;
0044         nWeighted += zWeights[ii - zProjections.begin()];
0045         centerW += (*ii) * zWeights[ii - zProjections.begin()];
0046       }
0047       centerW /= nWeighted;  //calculate the weighted peak center
0048       if (nWeighted > maxW) {
0049         maxW = nWeighted;
0050         centerWMax = centerW;
0051       }
0052     }
0053   }
0054   return centerWMax;
0055 }
0056
0057 #endif