Skip to content
Snippets Groups Projects
Commit 21bd1b26 authored by Chao Peng's avatar Chao Peng
Browse files

fix several typos and add more information in wfanalyzer test script

parent 34c196ac
No related branches found
No related tags found
No related merge requests found
......@@ -101,7 +101,7 @@ public:
WFData data;
if (!nsamples) { return data; }
auto buffer = ResSpectrum(samples, nsamples, res);
auto buffer = SmoothSpectrum(samples, nsamples, res);
// search local maxima
std::vector<Peak> candidates;
......@@ -171,30 +171,37 @@ public:
peak.integral = buffer[peak.pos] - data.ped.mean;
for (size_t i = peak.pos - 1; i >= peak.left; --i) {
double val = buffer[i] - data.ped.mean;
// stop when it touches the baseline
if (std::abs(val) < data.ped.err) { break; }
// stop when it touches or acrosses the baseline
if (std::abs(val) < data.ped.err || val * sample_height < 0.) {
peak.left = i; break;
}
peak.integral += val;
}
for (size_t i = peak.pos + 1; i <= peak.right; ++i) {
double val = buffer[i] - data.ped.mean;
if (std::abs(val) < data.ped.err) { break; }
if (std::abs(val) < data.ped.err || val * sample_height < 0.) {
peak.right = i; break;
}
peak.integral += val;
}
// determine the real sample peak
peak.height = samples[peak.pos];
peak.height = sample_height;
auto update_peak = [] (Peak &peak, double val, size_t pos) {
if (std::abs(val) > std::abs(peak.height)) {
peak.pos = pos;
peak.height = val;
}
};
size_t sample_pos = peak.pos;
for (size_t i = 1; i < res; ++i) {
if (samples[sample_pos - i] > peak.height) {
peak.pos = sample_pos - i;
peak.height = samples[sample_pos - i];
if (sample_pos > i) {
update_peak(peak, samples[sample_pos - i] - data.ped.mean, sample_pos - i);
}
if (samples[sample_pos + i] > peak.height) {
peak.pos = sample_pos + i;
peak.height = samples[sample_pos + i];
if (sample_pos + i < nsamples) {
update_peak(peak, samples[sample_pos + i] - data.ped.mean, sample_pos + i);
}
}
peak.height -= data.ped.mean;
// fill to results
data.peaks.emplace_back(peak);
}
......@@ -214,15 +221,15 @@ private:
public:
// static methods
template<typename T>
static std::vector<double> ResSpectrum(T *samples, size_t nsamples, size_t res)
static std::vector<double> SmoothSpectrum(T *samples, size_t nsamples, size_t res)
{
std::vector<double> buffer(nsamples);
for (size_t i = 0; i < nsamples; ++i) {
double val = samples[i];
double weights = 1.0;
for (size_t j = 0; j < res; ++j) {
for (size_t j = 1; j < res; ++j) {
if (j >= i || j + i >= nsamples) { continue; }
double weight = j/static_cast<double>(res + 1);
double weight = 1.0 - j/static_cast<double>(res + 1);
val += weight*(samples[i - j] + samples[i + j]);
weights += 2.*weight;
}
......
......@@ -7,8 +7,8 @@ void test(size_t res = 3, double thres = 10.0)
std::vector<std::vector<uint32_t>> tests = {
{
129, 132, 131, 128, 128, 127, 130, 130, 126, 126, 128, 126, 126,
127, 128, 125, 130, 130, 131, 143, 146, 149, 143, 142, 139, 136,
// 127, 128, 125, 123, 120, 121, 113, 109, 102, 109, 110, 119, 126,
// 127, 128, 125, 130, 130, 131, 143, 146, 149, 143, 142, 139, 136,
127, 128, 125, 123, 120, 121, 113, 109, 102, 109, 110, 119, 126,
134, 135, 132, 132, 132, 128, 139, 167, 188, 189, 180, 171, 164,
156, 150, 145, 143, 140, 136, 137, 136, 132, 131, 132, 132, 132,
130, 129, 131, 130, 131, 131, 128, 130, 129, 127, 128, 127
......@@ -29,18 +29,26 @@ void test(size_t res = 3, double thres = 10.0)
},
{
162, 158, 157, 159, 159, 160, 159, 161, 159, 160, 158, 159, 183,
// 236, 263, 271, 267, 254, 240, 228, 213, 205, 192, 186, 181, 177,
236, 1263, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 2271, 1177,
236, 263, 271, 267, 254, 240, 228, 213, 205, 192, 186, 181, 177,
// 236, 1263, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 8000, 2271, 1177,
174, 170, 165, 170, 165, 165, 164, 162, 161, 163, 163, 162, 159,
159, 161, 159, 161, 160, 158, 157, 159, 158, 156, 160, 159, 158,
160, 158, 158, 156, 158, 157, 159, 157, 159, 158, 156, 157
},
{
152, 150, 153, 151, 151, 154, 152, 151, 151, 151, 150, 154, 151,
154, 151, 155, 152, 152, 150, 152, 153, 153, 151, 150, 153, 151,
151, 150, 151, 153, 192, 244, 267, 266, 254, 238, 220, 206, 198,
188, 181, 174, 171, 167, 164, 161, 160, 160, 156, 157, 156, 156,
155, 155, 156, 154, 152, 153, 153, 190, 246, 269, 263, 250
},
};
wfa::Analyzer ana(res, thres);
auto c = new TCanvas("test", "test", 1280, 720);
c->DivideSquare(tests.size());
auto c = new TCanvas("test", "test", 1920, 1080);
c->DivideSquare(tests.size() + 1);
auto legend = new TLegend(0.1, 0.9, 0.9, 0.1);
int count = 0;
for (auto &samples : tests) {
......@@ -48,18 +56,19 @@ void test(size_t res = 3, double thres = 10.0)
c->cd(count);
// waveform
auto wf = new TGraph(samples.size());
wf->SetTitle("Waveform Analysis; Sample Number; ADC Values");
for (size_t i = 0; i < samples.size(); ++i) {
wf->SetPoint(i, i, samples[i]);
}
wf->SetLineColor(kBlue + 1);
wf->SetLineWidth(1);
wf->SetLineColor(kRed + 1);
wf->SetLineWidth(3);
wf->SetLineStyle(2);
wf->Draw("AL");
wf->GetXaxis()->SetRangeUser(0, samples.size() - 1);
// waveform resolution
auto wf2 = new TGraph(samples.size());
auto buf = wfa::Analyzer::ResSpectrum(&samples[0], samples.size(), res);
auto buf = wfa::Analyzer::SmoothSpectrum(&samples[0], samples.size(), res);
for (size_t i = 0; i < buf.size(); ++i) {
wf2->SetPoint(i, i, buf[i]);
}
......@@ -70,16 +79,38 @@ void test(size_t res = 3, double thres = 10.0)
// peak finding
auto data = ana.Analyze(&samples[0], samples.size());
auto ped = data.ped;
auto gr = new TGraph(samples.size());
gr->SetMarkerStyle(23);
for (auto &peak : data.peaks) {
gr->SetPoint(gr->GetN(), peak.pos, peak.height + data.ped.mean);
for (size_t i = 0; i < data.peaks.size(); ++i) {
auto color = i + 40;
auto peak = data.peaks[i];
auto grm = new TGraph(samples.size());
grm->SetMarkerStyle(peak.height > 0 ? 23 : 22);
grm->SetMarkerSize(1.5);
grm->SetMarkerColor(kBlack);
double range = wf->GetYaxis()-> GetXmax() - wf->GetYaxis()-> GetXmin();
double height = peak.height + data.ped.mean + (peak.height > 0 ? 1. : -1.5)*range*0.02;
grm->SetPoint(0, peak.pos, height);
grm->Draw("P same");
auto nint = peak.right - peak.left + 1;
auto grs = new TGraph(2*nint);
for (size_t i = 0; i < nint; ++i) {
auto val = buf[i + peak.left];
grs->SetPoint(i, i + peak.left, val);
grs->SetPoint(nint + i, peak.right - i, ped.mean);
}
grs->SetFillColor(color);
grs->SetFillStyle(3002);
grs->Draw("f same");
// 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.left << ", " << peak.right
<< ", " << peak.height << ", " << peak.integral << std::endl;
if (i == 0 && count == 1) {
legend->AddEntry(grm, Form("Peaks (threhold = %.2f)", thres), "p");
legend->AddEntry(grs, "Peak Integrals", "f");
}
}
gr->Draw("P same");
// pedestal line
/*
......@@ -101,11 +132,20 @@ void test(size_t res = 3, double thres = 10.0)
gre->SetPoint(i + samples.size(), j, pedy2 - ped.err);
}
gre->SetFillStyle(3001);
gre->SetFillColor(kRed);
gre->SetFillColor(kBlack);
grp->SetFillStyle(3001);
grp->SetFillColor(kBlack);
grp->SetLineStyle(0);
grp->SetLineWidth(2);
grp->SetLineColor(kRed);
grp->SetLineColor(kBlack);
gre->Draw("f same");
grp->Draw("L same");
if (count == 1) {
legend->AddEntry(wf, "Raw Samples", "l");
legend->AddEntry(wf2, Form("Smoothed Samples (res = %zu)", res), "l");
legend->AddEntry(grp, "Pedestal", "lf");
}
}
c->cd(tests.size() + 1);
legend->Draw();
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment