PeakFinder.cpp 4.2 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"

frey_m's avatar
frey_m committed
10 11 12 13 14 15 16 17 18 19 20
PeakFinder::PeakFinder(std::string elem, double min,
                       double max, double binWidth, bool singlemode)
    : element_m(elem)
    , binWidth_m(binWidth)
    , min_m(min)
    , max_m(max)
    , turn_m(0)
    , peakRadius_m(0.0)
    , registered_m(0)
    , singlemode_m(singlemode)
    , first_m(true)
21
    , finished_m(false)
22 23 24 25 26 27 28
{
    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
29 30


31
void PeakFinder::addParticle(const Vector_t& R) {
frey_m's avatar
frey_m committed
32
    
33
    double radius = std::hypot(R(0),R(1));
frey_m's avatar
frey_m committed
34
    radius_m.push_back(radius);
frey_m's avatar
frey_m committed
35
    
36 37 38 39 40 41 42
    peakRadius_m += radius;
    ++registered_m;
}


void PeakFinder::evaluate(const int& turn) {
    
frey_m's avatar
frey_m committed
43 44 45 46 47 48
    if ( first_m ) {
        turn_m = turn;
        first_m = false;
    }
    
    if ( turn_m != turn ) {
49
        finished_m = true;
frey_m's avatar
frey_m committed
50 51
    }
    
52 53
    bool globFinished = false;
    
54
    if ( !singlemode_m )
55 56 57
        allreduce(finished_m, globFinished, 1, std::logical_and<bool>());
    else
        globFinished = finished_m;
58
    
59
    if ( globFinished ) {
60
        
61 62
        this->computeCentroid_m();
        
63 64
        turn_m = turn;
        
65
        // reset
66 67
        peakRadius_m = 0.0;
        registered_m = 0;
68 69 70 71 72
        finished_m = false;
    }
}


frey_m's avatar
frey_m committed
73 74
void PeakFinder::save() {
    
75
    createHistogram_m();
76
    
frey_m's avatar
frey_m committed
77 78
    // last turn is not yet computed
    this->computeCentroid_m();
79
    
80
    if ( !peaks_m.empty() ) {
81 82
        // only rank 0 will go in here
        
83
        fn_m   = element_m + std::string(".peaks");
84 85 86 87 88 89 90 91 92 93 94 95 96 97
        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();
        
    }
98

frey_m's avatar
frey_m committed
99 100
    radius_m.clear();
    globHist_m.clear();
frey_m's avatar
frey_m committed
101 102
}

103

frey_m's avatar
frey_m committed
104 105 106
void PeakFinder::computeCentroid_m() {
    double globPeakRadius = 0.0;
    int globRegister = 0;
frey_m's avatar
frey_m committed
107
    
frey_m's avatar
frey_m committed
108 109
    //FIXME inefficient
    if ( !singlemode_m ) {
110 111
        reduce(peakRadius_m, globPeakRadius, 1, std::plus<double>());
        reduce(registered_m, globRegister, 1, std::plus<int>());
frey_m's avatar
frey_m committed
112
    } else {
113 114
        globPeakRadius = peakRadius_m;
        globRegister = registered_m;
frey_m's avatar
frey_m committed
115
    }
frey_m's avatar
frey_m committed
116
    
frey_m's avatar
frey_m committed
117
    if ( Ippl::myNode() == 0 ) {
118 119
        if ( globRegister > 0 )
            peaks_m.push_back(globPeakRadius / double(globRegister));
120
    }
frey_m's avatar
frey_m committed
121
}
frey_m's avatar
frey_m committed
122

123 124 125 126
void PeakFinder::createHistogram_m() {
    /*
     * create local histograms
     */
frey_m's avatar
frey_m committed
127
    
128 129
    globHist_m.resize(nBins_m);
    container_t locHist(nBins_m,0.0);
130

131 132
    double invBinWidth = 1.0 / binWidth_m;
    for(container_t::iterator it = radius_m.begin(); it != radius_m.end(); ++it) {
133 134
        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
135 136 137 138 139 140
        ++locHist[bin];
    }
    
    /*
     * create global histograms
     */
frey_m's avatar
frey_m committed
141 142 143 144 145 146 147 148
    if ( !singlemode_m ) {
        reduce(&(locHist[0]), &(locHist[0]) + locHist.size(),
               &(globHist_m[0]), OpAddAssign());
    } else {
        globHist_m.swap(locHist);
    }
    
//     reduce(locHist.data(), globHist_m.data(), locHist.size(), std::plus<double>());
149 150 151
}


frey_m's avatar
frey_m committed
152
void PeakFinder::open_m() {
153 154
    os_m.open(fn_m.c_str(), std::ios::out);
    hos_m.open(hist_m.c_str(), std::ios::out);
frey_m's avatar
frey_m committed
155 156 157 158
}


void PeakFinder::append_m() {
159 160
    os_m.open(fn_m.c_str(), std::ios::app);
    hos_m.open(hist_m.c_str(), std::ios::app);
frey_m's avatar
frey_m committed
161 162 163 164
}


void PeakFinder::close_m() {
165 166
    os_m.close();
    hos_m.close();
frey_m's avatar
frey_m committed
167 168 169 170
}


void PeakFinder::saveASCII_m() {
171 172 173
    os_m << "# Peak Radii (mm)" << std::endl;
    for (auto &radius : peaks_m)
        os_m << radius << std::endl;
frey_m's avatar
frey_m committed
174
        
175 176 177 178 179 180 181
    hos_m << "# Histogram bin counts (min, max, nbins, binsize) "
          << min_m << " mm "
          << max_m << " mm "
          << nBins_m << " "
          << binWidth_m << " mm" << std::endl;
    for (auto binCount : globHist_m)
        hos_m << binCount << std::endl;
frey_m's avatar
frey_m committed
182
}