PeakFinder.cpp 10.9 KB
Newer Older
frey_m's avatar
frey_m committed
1
#include "PeakFinder.h"
snuverink_j's avatar
snuverink_j committed
2 3 4

#include <algorithm>
#include <cmath>
frey_m's avatar
frey_m committed
5 6
#include <iterator>

frey_m's avatar
frey_m committed
7 8 9
#include "AbstractObjects/OpalData.h"
#include "Ippl.h"

10 11
PeakFinder::PeakFinder(std::string elem, double min, double max, double binWidth, bool singlemode):
    element_m(elem), binWidth_m(binWidth),
12 13 14 15 16 17 18 19
    min_m(min), max_m(max), singlemode_m(singlemode)
{
    if (min_m > max_m) {
        std::swap(min_m, max_m);
    }
    // calculate bins, round up so that histogram is large enough (add one for safety)
    nBins_m = static_cast<unsigned int>(std::ceil(( max_m - min_m ) / binWidth_m)) + 1;
}
frey_m's avatar
frey_m committed
20 21 22


void PeakFinder::addParticle(const Vector_t& R) {
frey_m's avatar
frey_m committed
23
    double radius = std::sqrt( dot(R, R) );
frey_m's avatar
frey_m committed
24 25 26 27 28 29
    radius_m.push_back(radius);
}


void PeakFinder::save() {
    
30
    createHistogram_m();
31

32 33
    bool found = false;
    
34
    if ( Ippl::getNodes() == 1 && singlemode_m == true) {
35 36 37 38 39 40 41
        found = findPeaks();
    } else {
        found = findPeaks(smoothingNumber_m,
                          minArea_m,
                          minFractionalArea_m,
                          minAreaAboveNoise_m,
                          minSlope_m);
42
    }
43
    
44
    if ( found ) {
45
        fn_m   = element_m + std::string(".peaks");
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
        hist_m = element_m + std::string(".hist");
        
        INFOMSG("Save " << fn_m << " and " << hist_m << endl);
        
        if(OpalData::getInstance()->inRestartRun())
            this->append_m();
        else
            this->open_m();
        
        this->saveASCII_m();
        
        this->close_m();
        Ippl::Comm->barrier();
        
    }
frey_m's avatar
frey_m committed
61
    
frey_m's avatar
frey_m committed
62 63
    radius_m.clear();
    globHist_m.clear();
frey_m's avatar
frey_m committed
64 65
    peakRadii_m.clear();
    fourSigmaPeaks_m.clear();
frey_m's avatar
frey_m committed
66 67
}

68

69 70 71 72 73 74 75
bool PeakFinder::findPeaks() {
    for (const auto &radius : radius_m)
        peakRadii_m.push_back(radius);
    return !peakRadii_m.empty();
}


76
bool PeakFinder::findPeaks(int smoothingNumber,
frey_m's avatar
frey_m committed
77 78 79 80
                           double minAreaFactor,
                           double minFractionalAreaFactor,
                           double minAreaAboveNoise,
                           double minSlope)
frey_m's avatar
frey_m committed
81
{
82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
    /* adapted from subroutine SEPAPR
     * Die Routine waehlt einen Beobachtungsindex.
     * Von diesem Aus wird fortlaufend die Peakflaeche
     * FTP integriert und mit dem aus dem letzten Messwert
     * und dessen Abstand vom Beobachtungspunkt gebildete
     * Dreieck ZPT verglichen. Ist FTP > ZPT ist ein neuer
     * Peak identifiziert. Der Beobachtungspunkt
     * verschiebt sich zum letzten Messwertindex und ab da
     * weiter, solange der Messwert abnimmt. Parallel wird
     * die Gesamtmessung aufintegriert, als 
     * zusaetzliches Kriterium zur Unterscheidung von echten
     * Peaks und Rauscheffekten.
     * 
     * smoothingNumber          Startindex in VAL für die Peakidentifikation
     * minAreaFactor            Zulässiger minimaler Anteil eines Einzelpeaks
     *                          am Messdatenintegral = Gewichtsfaktor für die
     *                          Elimination von Rauschpeaks
     * minFractionalAreaFactor  Gewichtsfaktor für die Gegenueberstellung FTP - ZPT
     *                          smoothen the data by summing neighbouring bins
     */
frey_m's avatar
frey_m committed
102
  
frey_m's avatar
frey_m committed
103
    const int size = static_cast<int>(globHist_m.size());
104
    if (size < smoothingNumber) {
105
        // no peaks can be found
106
        return false;
107
    }
frey_m's avatar
frey_m committed
108 109
    
    container_t smoothValues, sumSmoothValues;
110 111 112
    smoothValues.resize(size);
    sumSmoothValues.resize(size);
    double totalSum = 0.0;
frey_m's avatar
frey_m committed
113 114
    
    for (int i = smoothingNumber; i < size-smoothingNumber; i++) {
115
        double sum = 0.0;
frey_m's avatar
frey_m committed
116 117
        for (int j = -smoothingNumber; j <= smoothingNumber; j++) {
            sum += globHist_m[i+j];
118
        }
frey_m's avatar
frey_m committed
119
        sum /= smoothingNumber * 2 + 1;
120 121 122
        totalSum += sum;
        smoothValues[i] = sum;
        sumSmoothValues[i] = totalSum;
123
    }
frey_m's avatar
frey_m committed
124
    
125 126
    // set first and last values to value of smoothingNumber
    for (int i=0; i<smoothingNumber; i++) {
127 128
        smoothValues[i]        = smoothValues[smoothingNumber];
        smoothValues[size-1-i] = smoothValues[size-1-smoothingNumber];
frey_m's avatar
frey_m committed
129 130
    }
  
131 132
    std::vector<int> peakSeparatingIndices; // indices at minima (one more than number of peaks(!))
    peakSeparatingIndices.push_back(0);
frey_m's avatar
frey_m committed
133
  
134 135
    int nrPeaks            = 0;
    const double minArea   = minAreaFactor * totalSum; // minimum area for a peak
136
#ifdef PEAK_DEBUG
137
    INFOMSG("minArea " << minArea << endl);
138
#endif
139 140 141 142
    // number of indices corresponding to 10 mm
    const int maxIndex     = static_cast<int> (10 * size / binWidth_m);
    bool upwards           = false;
    bool newPeak           = false;
frey_m's avatar
frey_m committed
143 144
    
    for (int i = 1; i < size; i++) {
145 146 147 148 149
        int startIndex = std::max(i-maxIndex, peakSeparatingIndices.back());
        double ftp     = sumSmoothValues[i] - sumSmoothValues[startIndex];
        double ftpPeak = ftp - (i - startIndex)*smoothValues[startIndex]; // peak - noiselevel
        double slope   = (smoothValues[i] - smoothValues[startIndex]) / (i-startIndex);
        double zpt     = minFractionalAreaFactor * (smoothValues[i] - smoothValues[startIndex]) * (i - startIndex);
frey_m's avatar
frey_m committed
150
        
151
        if (ftpPeak >= zpt && ftp > minArea && ftpPeak > minAreaAboveNoise && slope > minSlope) {
152
#ifdef PEAK_DEBUG
153 154 155 156 157 158 159
            if (newPeak == false) {
                INFOMSG("Peak "     << peakSeparatingIndices.size() << endl);
                INFOMSG("Fraction " << ftpPeak << " " << zpt << endl);
                INFOMSG("Area "     << ftp     << " " << minArea << endl);
                INFOMSG("Noise "    << ftpPeak << " " << minAreaAboveNoise << endl);
                INFOMSG("Slope "    << slope   << " " << minSlope << endl);
            }
160
#endif
161 162
            newPeak = true;
        }
frey_m's avatar
frey_m committed
163
        
164 165 166 167 168 169 170 171 172 173 174 175 176 177
        if (smoothValues[i] > smoothValues [i-1] || i == size-1) {
            if (upwards == false || i == size-1) {
                upwards = true;
                if (newPeak == true) {
                    nrPeaks++;
                    peakSeparatingIndices.push_back(i-1);
                    newPeak = false;
                } else if (smoothValues[peakSeparatingIndices.back()] >= smoothValues[i]) {
                    peakSeparatingIndices.back() = i;
                }
            }
        } else {
            upwards = false;
        }
178
    }
frey_m's avatar
frey_m committed
179
    
180
    // debug
181
#ifdef PEAK_DEBUG
182
    INFOMSG("Number of peaks found: " << nrPeaks << endl);
183
#endif
184 185
    peakRadii_m.resize(nrPeaks);
    fourSigmaPeaks_m.resize(nrPeaks);
186
    
frey_m's avatar
frey_m committed
187
    // the position of the peak is bin centered
188 189
    container_t positions;
    positions.reserve(nBins_m);
frey_m's avatar
frey_m committed
190 191
    
    for (unsigned int i = 0; i < nBins_m; i++) {
192
        positions.push_back(min_m + (i + 0.5) * binWidth_m);
193 194
    }
    
frey_m's avatar
frey_m committed
195
    for (int i = 1; i < (int)(peakSeparatingIndices.size()); i++) {
196 197
        int startIndex = peakSeparatingIndices[i-1];
        int endIndex   = peakSeparatingIndices[i];
frey_m's avatar
frey_m committed
198 199
        analysePeak(globHist_m, positions, startIndex, endIndex,
                    peakRadii_m[i-1], fourSigmaPeaks_m[i-1]);
frey_m's avatar
frey_m committed
200
    }
201 202
    
    return !peakRadii_m.empty();
frey_m's avatar
frey_m committed
203 204
}

205

frey_m's avatar
frey_m committed
206
void PeakFinder::analysePeak(const container_t& values,
207 208 209 210
                             const container_t& positions,
                             const int startIndex, const int endIndex,
                             double& peak,
                             double& fourSigma)const
frey_m's avatar
frey_m committed
211
{
212 213 214 215 216 217
    // original subroutine ANALPR
    int range      = endIndex - startIndex;
    // find maximum
    double maximum   = -1;
    int maximumIndex = -1;
    int relMaxIndex  = -1;
frey_m's avatar
frey_m committed
218 219
    
    for (int j = startIndex; j <= endIndex; j++) {
220 221 222 223 224
        if (values[j] > maximum) {
            maximum = values[j];
            maximumIndex = j;
            relMaxIndex  = j - startIndex; // count from peak separation
        }
frey_m's avatar
frey_m committed
225
    }
frey_m's avatar
frey_m committed
226
    
227
    peak = positions[maximumIndex];
frey_m's avatar
frey_m committed
228
    
229 230 231
    // left limits, go down from peak to separation index
    int index20 = -1;
    int indexLeftEnd = 0; // left limit of peak
frey_m's avatar
frey_m committed
232 233
    
    for (int j = relMaxIndex; j >= 0; j--) {
234 235
        int index = j + startIndex;
        double value = values[index];
frey_m's avatar
frey_m committed
236 237 238
        if (value > 0.2 * maximum)
            index20 = j; // original code had i-1
        
239
        // if too far out, then break (not sure where formula comes from)
frey_m's avatar
frey_m committed
240 241
        if ( j < (3 * index20 - 2 * relMaxIndex) ) {
            indexLeftEnd = j;
242 243
            break;
        }
frey_m's avatar
frey_m committed
244
    }
frey_m's avatar
frey_m committed
245
    
246 247 248
    // right limits
    index20 = -1;
    int indexRightEnd = range; // right limit of peak
frey_m's avatar
frey_m committed
249
    
250
    // loop on right side of peak
frey_m's avatar
frey_m committed
251
    for (int j = relMaxIndex; j <= range; j++) {
252 253 254 255
        int index = j + startIndex;
        double value = values[index];
        if (value > 0.2 *maximum) {index20    = j;}
        // if too far out, then break (not sure where formula comes from)
frey_m's avatar
frey_m committed
256 257
        if ( j > (3 * index20 - 2 * relMaxIndex) ) {
            indexRightEnd = j;
258 259
            break;
        }
frey_m's avatar
frey_m committed
260 261
    }
    
262
    if (indexRightEnd - indexLeftEnd == 0) { // no peak
263 264
        fourSigma = 0.0;
        return; // return zeros for sigma
265
    }
frey_m's avatar
frey_m committed
266
    
267
    double sum=0.0, radialSum=0.0;
frey_m's avatar
frey_m committed
268 269
    
    for (int j = indexLeftEnd; j <= indexRightEnd; j++) {
270 271 272
        int index = j + startIndex;
        sum       += values[index];
        radialSum += values[index] * positions[index];
273
    }
frey_m's avatar
frey_m committed
274
    
275 276
    double mean = radialSum / sum;
    double variance = 0.0;
frey_m's avatar
frey_m committed
277 278
    
    for (int j = indexLeftEnd; j <= indexRightEnd; j++) {
279 280 281 282
        int index = j + startIndex;
        double value = values[index];
        double dx = positions[index] - mean;
        variance += value * dx * dx;
283 284
    }
    fourSigma = 4 * std::sqrt(variance / sum);
frey_m's avatar
frey_m committed
285
}
frey_m's avatar
frey_m committed
286 287


288 289 290 291
void PeakFinder::createHistogram_m() {
    /*
     * create local histograms
     */
frey_m's avatar
frey_m committed
292
    
293 294
    globHist_m.resize(nBins_m);
    container_t locHist(nBins_m,0.0);
295

296 297
    double invBinWidth = 1.0 / binWidth_m;
    for(container_t::iterator it = radius_m.begin(); it != radius_m.end(); ++it) {
298 299
        int bin = static_cast<int>(std::abs(*it - min_m ) * invBinWidth);
        if (bin < 0 || (unsigned int)bin >= nBins_m) continue; // Probe might save particles outside its boundary
300 301 302 303 304 305
        ++locHist[bin];
    }
    
    /*
     * create global histograms
     */
frey_m's avatar
frey_m committed
306 307
    reduce(&(locHist[0]), &(locHist[0]) + locHist.size(),
           &(globHist_m[0]), OpAddAssign());
308 309 310
}


frey_m's avatar
frey_m committed
311 312 313
void PeakFinder::open_m() {
    if ( Ippl::myNode() == 0 ) {
        os_m.open(fn_m.c_str(), std::ios::out);
frey_m's avatar
frey_m committed
314
        hos_m.open(hist_m.c_str(), std::ios::out);
frey_m's avatar
frey_m committed
315 316 317 318 319 320 321
    }
}


void PeakFinder::append_m() {
    if ( Ippl::myNode() == 0 ) {
        os_m.open(fn_m.c_str(), std::ios::app);
frey_m's avatar
frey_m committed
322
        hos_m.open(hist_m.c_str(), std::ios::app);
frey_m's avatar
frey_m committed
323 324 325 326 327
    }
}


void PeakFinder::close_m() {
frey_m's avatar
frey_m committed
328
    if ( Ippl::myNode() == 0 ) {
frey_m's avatar
frey_m committed
329
        os_m.close();
frey_m's avatar
frey_m committed
330 331
        hos_m.close();
    }
frey_m's avatar
frey_m committed
332 333 334 335
}


void PeakFinder::saveASCII_m() {
frey_m's avatar
frey_m committed
336
    if ( Ippl::myNode() == 0 )  {
frey_m's avatar
frey_m committed
337 338
        os_m << "# Peak Radii (mm)" << std::endl;
        for (auto &radius : peakRadii_m)
339
            os_m << radius << std::endl;
frey_m's avatar
frey_m committed
340
        
341
        hos_m << "# Histogram bin counts (min, max, nbins, binsize) "
342 343
              << min_m << " mm "
              << max_m << " mm "
344
              << nBins_m << " "
345
              << binWidth_m << " mm" << std::endl;
frey_m's avatar
frey_m committed
346 347
        for (auto binCount : globHist_m)
            hos_m << binCount << std::endl;
frey_m's avatar
frey_m committed
348 349 350
    }
}