Commit b8fe4db5 authored by snuverink_j's avatar snuverink_j
Browse files

peak finder histogram gets fixed size (of probe), add singlemode boolean at initialisation

parent eaf072ea
......@@ -63,7 +63,7 @@ public:
FORCE = 1,
AUTO = 2
};
typedef std::vector<double> dvector_t;
typedef std::vector<int> ivector_t;
......
......@@ -99,17 +99,23 @@ void Probe::initialise(PartBunchBase<double, 3> *bunch, double &startField, doub
void Probe::initialise(PartBunchBase<double, 3> *bunch) {
*gmsg << "* Initialize probe" << endl;
if (filename_m == std::string("")) {
peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder(getName()));
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
} else {
peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder(filename_m.substr(0, filename_m.rfind("."))));
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
}
double r_start = std::sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = std::sqrt(xend_m * xend_m + yend_m * yend_m);
std::string name;
if (filename_m == std::string(""))
name = getName();
else
name = filename_m.substr(0, filename_m.rfind("."));
bool singlemode = (bunch->getTotalNum() == 1) ? true : false;
peakfinder_m = std::unique_ptr<PeakFinder> (new PeakFinder (name, r_start, r_end, singlemode));
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(name, !Options::asciidump));
}
void Probe::finalise() {
*gmsg << "* Finalize probe " << getName() << endl;
*gmsg << "* Finalize probe " << getName() << endl;
peakfinder_m->save();
lossDs_m->save();
}
......
......@@ -7,16 +7,16 @@
#include "AbstractObjects/OpalData.h"
#include "Ippl.h"
#include "Algorithms/PartBunchBase.h"
PeakFinder::PeakFinder(std::string elem):
element_m(elem), nBins_m(10), binWidth_m(1.0 /*mm*/),
globMin_m(0.0),globMax_m(5000.0)
{}
PeakFinder::PeakFinder() : PeakFinder(std::string("NULL"))
{ }
PeakFinder::PeakFinder(std::string elem, double min, double max, bool singlemode):
element_m(elem), binWidth_m(1.0 /*mm*/),
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;
}
void PeakFinder::addParticle(const Vector_t& R) {
......@@ -31,9 +31,7 @@ void PeakFinder::save() {
bool found = false;
std::size_t nParticles = (OpalData::getInstance()->getPartBunch())->getTotalNum();
if ( Ippl::getNodes() == 1 && nParticles == 1) {
if ( Ippl::getNodes() == 1 && singlemode_m == true) {
found = findPeaks();
} else {
found = findPeaks(smoothingNumber_m,
......@@ -41,7 +39,7 @@ void PeakFinder::save() {
minFractionalArea_m,
minAreaAboveNoise_m,
minSlope_m);
}
}
if ( found ) {
fn_m = element_m + std::string(".peaks");
......@@ -191,7 +189,7 @@ bool PeakFinder::findPeaks(int smoothingNumber,
positions.reserve(nBins_m);
for (unsigned int i = 0; i < nBins_m; i++) {
positions.push_back(globMin_m + (i + 0.5) * binWidth_m);
positions.push_back(min_m + (i + 0.5) * binWidth_m);
}
for (int i = 1; i < (int)(peakSeparatingIndices.size()); i++) {
......@@ -288,42 +286,17 @@ void PeakFinder::analysePeak(const container_t& values,
void PeakFinder::createHistogram_m() {
/* A core might have no particles, thus, using initial high values
* in order to find real global minimum and maximum among all cores.
* It might happen that no particles hit the probe and subsequently the
* container radius_m would be empty for all MPI processes. In that case we
* fix the number of bins to a small value in order to avoid a drastic memory
* increase.
*/
double locMin = 1e10, locMax = -1e10;
if (!radius_m.empty()) {
// compute global minimum and maximum radius
auto result = std::minmax_element(radius_m.begin(), radius_m.end());
locMin = *result.first;
locMax = *result.second;
}
reduce(locMin, globMin_m, OpMinAssign());
reduce(locMax, globMax_m, OpMaxAssign());
/*
* create local histograms
*/
if (globMax_m < -1e9)
nBins_m = 10; // no particles in probe
else {
// calculate bins, round up so that histogram is large enough (add one for safety)
nBins_m = static_cast<unsigned int>(std::ceil(( globMax_m - globMin_m ) / binWidth_m)) + 1;
}
globHist_m.resize(nBins_m);
container_t locHist(nBins_m,0.0);
double invBinWidth = 1.0 / binWidth_m;
for(container_t::iterator it = radius_m.begin(); it != radius_m.end(); ++it) {
int bin = static_cast<int>(std::abs(*it - globMin_m ) * invBinWidth);
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
++locHist[bin];
}
......@@ -366,8 +339,8 @@ void PeakFinder::saveASCII_m() {
os_m << radius << std::endl;
hos_m << "# Histogram bin counts (min, max, nbins, binsize) "
<< globMin_m << " mm "
<< globMax_m << " mm "
<< min_m << " mm "
<< max_m << " mm "
<< nBins_m << " "
<< binWidth_m << " mm" << std::endl;
for (auto binCount : globHist_m)
......
......@@ -28,9 +28,9 @@ public:
public:
PeakFinder();
PeakFinder() = delete;
explicit PeakFinder(std::string elem);
PeakFinder(std::string elem, double min, double max, bool singlemode);
/*!
* Append the particle coordinates to the container
......@@ -92,7 +92,7 @@ private:
const int startIndex, const int endIndex,
double& peak,
double& fourSigma)const;
private:
container_t radius_m;
/// global histogram values
......@@ -119,8 +119,10 @@ private:
/// Bin width in mm
double binWidth_m;
///@{ histogram size
double globMin_m, globMax_m;
double min_m, max_m;
///@}
/// Single particle mode on (different peakfinder)
bool singlemode_m;
///@{ Peak analysis parameters (copied from RRI2 probe program for now, need to be tuned a bit)
const int smoothingNumber_m = 0;
const double minArea_m = 0.025;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment