From 8046936493e952b75a479e62709be1ee90087b3e Mon Sep 17 00:00:00 2001
From: Chao Peng <cpeng@anl.gov>
Date: Tue, 4 Aug 2020 15:24:45 -0500
Subject: [PATCH] update wfanalyzer
---
include/WfAnalyzer.h | 49 ++++++++++++++++++++-------------------
tools/test_wfanalyzer.cxx | 9 +++----
2 files changed, 30 insertions(+), 28 deletions(-)
diff --git a/include/WfAnalyzer.h b/include/WfAnalyzer.h
index c6f5b7f..2e3a9fa 100644
--- a/include/WfAnalyzer.h
+++ b/include/WfAnalyzer.h
@@ -14,10 +14,6 @@
namespace wfa {
// data structures
-struct Pedestal {
- double mean, err;
-};
-
struct Peak {
size_t pos, left, right;
double height, integral;
@@ -27,6 +23,15 @@ struct Peak {
}
};
+struct Pedestal {
+ double mean, err;
+};
+
+struct WFData {
+ Pedestal ped;
+ std::vector<Peak> peaks;
+};
+
// calculate mean and standard deviation of an array
template<typename T>
inline void _calc_mean_err(double &mean, double &err, T *source, size_t npts)
@@ -91,10 +96,10 @@ public:
: thres(threshold), res(resolution) {}
~Analyzer() {}
- void Analyze(uint32_t *samples, size_t nsamples)
+ WFData Analyze(uint32_t *samples, size_t nsamples)
{
- peaks.clear();
- if (!nsamples) { return; }
+ WFData data;
+ if (!nsamples) { return data; }
auto buffer = ResSpectrum(samples, nsamples, res);
@@ -109,7 +114,7 @@ public:
int tr1 = trend(buffer[i], buffer[i - 1]);
int tr2 = trend(buffer[i], buffer[i + 1]);
// peak at the rising (declining) edge
- if ((tr1 * tr2 >= 0) && (std::abs(tr1) > 0) && std::abs(buffer[i]) > ped.err*3.0) {
+ if ((tr1 * tr2 >= 0) && (std::abs(tr1) > 0)) {
size_t left = 1, right = 1;
// search the peak range
while ((i > left + 1) && (trend(buffer[i - left], buffer[i - left - 1]) == tr1)) {
@@ -147,51 +152,47 @@ public:
// calculate pedestal
if (nped <= nsamples/4) {
- ped = CalcPedestal(&buffer[0], nsamples, 1.5, 3);
+ data.ped = CalcPedestal(&buffer[0], nsamples, 1.5, 3);
} else {
- ped = CalcPedestal(&ybuf[0], nped, 1.2, 3);
+ data.ped = CalcPedestal(&ybuf[0], nped, 1.2, 3);
}
for (auto &peak : candidates) {
// real sample height
- double sample_height = samples[peak.pos] - ped.mean;
+ double sample_height = samples[peak.pos] - data.ped.mean;
if ((sample_height * peak.height < 0.) ||
(std::abs(sample_height) < thres) ||
- (std::abs(sample_height) < 3.0*ped.err)) {
+ (std::abs(sample_height) < 3.0*data.ped.err)) {
continue;
}
// integrate it over the peak range
peak.height = sample_height;
- peak.integral = buffer[peak.pos] - ped.mean;
+ peak.integral = buffer[peak.pos] - data.ped.mean;
for (size_t i = 1; i <= peak.left; ++i) {
- double val = buffer[peak.pos - 1] - ped.mean;
- if (std::abs(val) < ped.err) { break; }
+ double val = buffer[peak.pos - 1] - data.ped.mean;
+ // stop when it touches the baseline
+ if (std::abs(val) < data.ped.err) { break; }
peak.integral += val;
}
for (size_t i = 1; i <= peak.right; ++i) {
- double val = buffer[peak.pos + i] - ped.mean;
- if (std::abs(val) < ped.err) { break; }
+ double val = buffer[peak.pos + i] - data.ped.mean;
+ if (std::abs(val) < data.ped.err) { break; }
peak.integral += val;
}
- peaks.emplace_back(peak);
+ data.peaks.emplace_back(peak);
}
/*
std::sort(peaks.begin(), peaks.end(),
[] (const Peak &p1, const Peak &p2) { return p1.height > p2.height; });
*/
-
+ return data;
}
- const Pedestal &GetPedestal() const { return ped; }
- const std::vector<Peak> &GetPeaks() const { return peaks; }
-
private:
double thres;
size_t res;
- Pedestal ped;
- std::vector<Peak> peaks;
public:
diff --git a/tools/test_wfanalyzer.cxx b/tools/test_wfanalyzer.cxx
index 046ee82..e4c7fc9 100644
--- a/tools/test_wfanalyzer.cxx
+++ b/tools/test_wfanalyzer.cxx
@@ -40,6 +40,7 @@ void test(size_t res = 3, double thres = 10.0)
wf->SetLineWidth(1);
wf->SetLineStyle(2);
wf->Draw("AL");
+ wf->GetXaxis()->SetRangeUser(0, samples.size() - 1);
// waveform resolution
auto wf2 = new TGraph(samples.size());
@@ -52,12 +53,12 @@ void test(size_t res = 3, double thres = 10.0)
wf2->Draw("L same");
// peak finding
- ana.Analyze(&samples[0], samples.size());
- auto ped = ana.GetPedestal();
+ auto data = ana.Analyze(&samples[0], samples.size());
+ auto ped = data.ped;
auto gr = new TGraph(samples.size());
gr->SetMarkerStyle(23);
- for (auto &peak : ana.GetPeaks()) {
- gr->SetPoint(gr->GetN(), peak.pos, peak.height + ped.mean);
+ for (auto &peak : data.peaks) {
+ gr->SetPoint(gr->GetN(), peak.pos, peak.height + data.ped.mean);
// gr->SetPoint(gr->GetN(), peak.pos - peak.left, peak.height + peak.base);
// gr->SetPoint(gr->GetN(), peak.pos + peak.right, peak.height + peak.base);
std::cout << peak.pos << ", " << peak.pos - peak.left << ", " << peak.pos + peak.right
--
GitLab