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

improve negative peak finding and graph formats

parent 8969c9c6
Branches master
No related tags found
No related merge requests found
...@@ -111,13 +111,15 @@ public: ...@@ -111,13 +111,15 @@ public:
// get final results // get final results
for (auto &peak : candidates) { for (auto &peak : candidates) {
// real sample height // pedestal subtraction
double sample_height = samples[peak.pos] - data.ped.mean; double peak_height = buffer[peak.pos] - data.ped.mean;
if ((sample_height * peak.height < 0.) || // wrong baselin in the rough candidtes finding, below threshold, or not statistically significant
(std::abs(sample_height) < _thres) || if ((peak_height * peak.height < 0.) ||
(std::abs(sample_height) < 3.0*data.ped.err)) { (std::abs(peak_height) < _thres) ||
(std::abs(peak_height) < 3.0*data.ped.err)) {
continue; continue;
} }
peak.height = peak_height;
// integrate it over the peak range // integrate it over the peak range
peak.integral = buffer[peak.pos] - data.ped.mean; peak.integral = buffer[peak.pos] - data.ped.mean;
...@@ -125,28 +127,28 @@ public: ...@@ -125,28 +127,28 @@ public:
for (int i = peak.pos - 1; i >= static_cast<int>(peak.left); --i) { for (int i = peak.pos - 1; i >= static_cast<int>(peak.left); --i) {
double val = buffer[i] - data.ped.mean; double val = buffer[i] - data.ped.mean;
// stop when it touches or acrosses the baseline // stop when it touches or acrosses the baseline
if (std::abs(val) < data.ped.err || val * sample_height < 0.) { if (std::abs(val) < data.ped.err || val * peak.height < 0.) {
peak.left = i; break; peak.left = i; break;
} }
peak.integral += val; peak.integral += val;
} }
for (size_t i = peak.pos + 1; i <= peak.right; ++i) { for (size_t i = peak.pos + 1; i <= peak.right; ++i) {
double val = buffer[i] - data.ped.mean; double val = buffer[i] - data.ped.mean;
if (std::abs(val) < data.ped.err || val * sample_height < 0.) { if (std::abs(val) < data.ped.err || val * peak.height < 0.) {
peak.right = i; break; peak.right = i; break;
} }
peak.integral += val; peak.integral += val;
} }
// determine the real sample peak // determine the real sample peak
peak.height = sample_height; size_t sample_pos = peak.pos;
peak.height = samples[sample_pos] - data.ped.mean;
auto update_peak = [] (Peak &peak, double val, size_t pos) { auto update_peak = [] (Peak &peak, double val, size_t pos) {
if (std::abs(val) > std::abs(peak.height)) { if (std::abs(val) > std::abs(peak.height)) {
peak.pos = pos; peak.pos = pos;
peak.height = val; peak.height = val;
} }
}; };
size_t sample_pos = peak.pos;
for (size_t i = 1; i < _res; ++i) { for (size_t i = 1; i < _res; ++i) {
if (sample_pos > i) { if (sample_pos > i) {
update_peak(peak, samples[sample_pos - i] - data.ped.mean, sample_pos - i); update_peak(peak, samples[sample_pos - i] - data.ped.mean, sample_pos - i);
...@@ -189,7 +191,7 @@ public: ...@@ -189,7 +191,7 @@ public:
// progressively find a good baseline // progressively find a good baseline
ped.mean = max_mean; ped.mean = max_mean;
for (size_t i = 0; i <= buffer.size() - ntrails; ++i) { for (size_t i = 0; i <= buffer.size() - ntrails; ++i) {
double mean, err; double mean = 0., err = 100.*good_err;
_calc_mean_err(mean, err, &buffer[i], ntrails); _calc_mean_err(mean, err, &buffer[i], ntrails);
if(err < good_err && mean < max_mean) { if(err < good_err && mean < max_mean) {
find_baseline = true; find_baseline = true;
...@@ -272,7 +274,7 @@ public: ...@@ -272,7 +274,7 @@ public:
std::vector<Peak> candidates; std::vector<Peak> candidates;
if (buffer.size() < 3) { return candidates; } if (buffer.size() < 3) { return candidates; }
candidates.reserve(buffer.size()/2); candidates.reserve(buffer.size()/3);
// get trend // get trend
auto trend = [] (double v1, double v2, double thr = 0.1) { auto trend = [] (double v1, double v2, double thr = 0.1) {
return std::abs(v1 - v2) < thr ? 0 : (v1 > v2 ? 1 : -1); return std::abs(v1 - v2) < thr ? 0 : (v1 > v2 ? 1 : -1);
...@@ -291,8 +293,7 @@ public: ...@@ -291,8 +293,7 @@ public:
right ++; right ++;
} }
double base = buffer[i - left] * right / static_cast<double>(left + right) double base = (buffer[i - left] * right + buffer[i + right] * left) / static_cast<double>(left + right);
+ buffer[i + right] * left / static_cast<double>(left + right);
double height = std::abs(buffer[i] - base); double height = std::abs(buffer[i] - base);
if (height > height_thres) { if (height > height_thres) {
candidates.push_back(Peak{i, i - left, i + right, buffer[i] - base, 0}); candidates.push_back(Peak{i, i - left, i + right, buffer[i] - base, 0});
......
...@@ -28,7 +28,8 @@ struct WFGraph ...@@ -28,7 +28,8 @@ struct WFGraph
WFData result; WFData result;
}; };
WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &samples, bool show_tspec = false) WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &samples, bool show_tspec = false,
double shade_alpha = 0.3)
{ {
WFGraph res; WFGraph res;
...@@ -47,7 +48,7 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam ...@@ -47,7 +48,7 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam
wf->SetPoint(i, i, samples[i]); wf->SetPoint(i, i, samples[i]);
} }
wf->SetLineColor(kRed + 1); wf->SetLineColor(kRed + 1);
wf->SetLineWidth(3); wf->SetLineWidth(1);
wf->SetLineStyle(2); wf->SetLineStyle(2);
res.mg->Add(wf, "l"); res.mg->Add(wf, "l");
res.objs.push_back(dynamic_cast<TObject*>(wf)); res.objs.push_back(dynamic_cast<TObject*>(wf));
...@@ -60,7 +61,7 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam ...@@ -60,7 +61,7 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam
wf2->SetPoint(i, i, buf[i]); wf2->SetPoint(i, i, buf[i]);
} }
wf2->SetLineColor(kBlue + 1); wf2->SetLineColor(kBlue + 1);
wf2->SetLineWidth(3); wf2->SetLineWidth(2);
res.mg->Add(wf2, "l"); res.mg->Add(wf2, "l");
res.objs.push_back(dynamic_cast<TObject*>(wf2)); res.objs.push_back(dynamic_cast<TObject*>(wf2));
res.entries.emplace_back(LegendEntry{wf2, Form("Smoothed Samples (res = %zu)", ana.GetResolution()), "l"}); res.entries.emplace_back(LegendEntry{wf2, Form("Smoothed Samples (res = %zu)", ana.GetResolution()), "l"});
...@@ -92,8 +93,8 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam ...@@ -92,8 +93,8 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam
grs->SetPoint(i, i + peak.left, val); grs->SetPoint(i, i + peak.left, val);
grs->SetPoint(nint + i, peak.right - i, ped.mean); grs->SetPoint(nint + i, peak.right - i, ped.mean);
} }
grs->SetFillColor(color); grs->SetFillColorAlpha(color, shade_alpha);
grs->SetFillStyle(3001); grs->SetFillStyle(3000);
res.mg->Add(grs, "f"); res.mg->Add(grs, "f");
res.objs.push_back(dynamic_cast<TObject*>(grs)); res.objs.push_back(dynamic_cast<TObject*>(grs));
if (i == 0) { if (i == 0) {
...@@ -101,13 +102,13 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam ...@@ -101,13 +102,13 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam
} }
} }
grm_pos->SetMarkerStyle(23); grm_pos->SetMarkerStyle(23);
grm_pos->SetMarkerSize(1.5); grm_pos->SetMarkerSize(1.0);
if (grm_pos->GetN()) { res.mg->Add(grm_pos, "p"); } if (grm_pos->GetN()) { res.mg->Add(grm_pos, "p"); }
res.objs.push_back(dynamic_cast<TObject*>(grm_pos)); res.objs.push_back(dynamic_cast<TObject*>(grm_pos));
res.entries.emplace_back(LegendEntry{grm_pos, "Peaks", "p"}); res.entries.emplace_back(LegendEntry{grm_pos, "Peaks", "p"});
grm_neg->SetMarkerStyle(22); grm_neg->SetMarkerStyle(22);
grm_neg->SetMarkerSize(1.5); grm_neg->SetMarkerSize(1.0);
if (grm_neg->GetN()) { res.mg->Add(grm_neg, "p"); } if (grm_neg->GetN()) { res.mg->Add(grm_neg, "p"); }
res.objs.push_back(dynamic_cast<TObject*>(grm_neg)); res.objs.push_back(dynamic_cast<TObject*>(grm_neg));
...@@ -117,9 +118,8 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam ...@@ -117,9 +118,8 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam
grp->SetPoint(i, i, ped.mean); grp->SetPoint(i, i, ped.mean);
grp->SetPointError(i, 0, ped.err); grp->SetPointError(i, 0, ped.err);
} }
grp->SetFillColor(kBlack); grp->SetFillColorAlpha(kBlack, shade_alpha);
grp->SetFillStyle(3001); grp->SetFillStyle(3000);
grp->SetLineStyle(0);
grp->SetLineWidth(2); grp->SetLineWidth(2);
grp->SetLineColor(kBlack); grp->SetLineColor(kBlack);
res.mg->Add(grp, "l3"); res.mg->Add(grp, "l3");
......
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
// channel look-up table // channel look-up table
static const std::vector<std::vector<std::string>> _channels = { static const std::vector<std::vector<std::string>> _quads = {
{"Cer11_1", "Cer11_2", "Cer11_3", "Cer11_4"}, {"Cer11_1", "Cer11_2", "Cer11_3", "Cer11_4"},
{"Cer12_1", "Cer12_2", "Cer12_3", "Cer12_4"}, {"Cer12_1", "Cer12_2", "Cer12_3", "Cer12_4"},
{"Cer13_1", "Cer13_2", "Cer13_3", "Cer13_4"}, {"Cer13_1", "Cer13_2", "Cer13_3", "Cer13_4"},
...@@ -28,14 +28,15 @@ static const std::vector<std::vector<std::string>> _channels = { ...@@ -28,14 +28,15 @@ static const std::vector<std::vector<std::string>> _channels = {
{"Cer42_1", "Cer42_2", "Cer42_3", "Cer42_4"}, {"Cer42_1", "Cer42_2", "Cer42_3", "Cer42_4"},
{"Cer43_1", "Cer43_2", "Cer43_3", "Cer43_4"}, {"Cer43_1", "Cer43_2", "Cer43_3", "Cer43_4"},
{"Cer44_1", "Cer44_2", "Cer44_3", "Cer44_4"}, {"Cer44_1", "Cer44_2", "Cer44_3", "Cer44_4"},
/*
{"Cer11_5"}, {"Cer12_5"}, {"Cer13_5"}, {"Cer14_5"},
{"Cer21_5"}, {"Cer22_5"}, {"Cer23_5"}, {"Cer24_5"},
{"Cer31_5"}, {"Cer32_5"}, {"Cer33_5"}, {"Cer34_5"},
{"Cer41_5"}, {"Cer42_5"}, {"Cer43_5"}, {"Cer44_5"},
*/
}; };
static const std::vector<std::vector<std::string>> _sums = {{
"Cer11_5", "Cer12_5", "Cer13_5", "Cer14_5",
"Cer21_5", "Cer22_5", "Cer23_5", "Cer24_5",
"Cer31_5", "Cer32_5", "Cer33_5", "Cer34_5",
"Cer41_5", "Cer42_5", "Cer43_5", "Cer44_5",
}};
struct TestChannel { struct TestChannel {
std::string name; std::string name;
std::vector<uint32_t> raw; std::vector<uint32_t> raw;
...@@ -47,7 +48,8 @@ struct TestEvent { ...@@ -47,7 +48,8 @@ struct TestEvent {
}; };
// get events from root file // get events from root file
std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indices = {}, int nev = 5) std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indices = {}, int nev = 5,
const std::vector<std::vector<std::string>> &channels = _quads)
{ {
std::vector<TestEvent> res; std::vector<TestEvent> res;
// root tree // root tree
...@@ -56,13 +58,13 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi ...@@ -56,13 +58,13 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi
int maxn = t->GetEntries(); int maxn = t->GetEntries();
// fetch tree data // fetch tree data
std::vector<std::vector<int>> nraws(_channels.size()); std::vector<std::vector<int>> nraws(channels.size());
std::vector<std::vector<std::vector<int>>> buffers(_channels.size()); std::vector<std::vector<std::vector<int>>> buffers(channels.size());
for (size_t i = 0; i < _channels.size(); ++i) { for (size_t i = 0; i < channels.size(); ++i) {
nraws[i].resize(_channels[i].size()); nraws[i].resize(channels[i].size());
buffers[i].resize(_channels[i].size()); buffers[i].resize(channels[i].size());
for (size_t j = 0; j < _channels[i].size(); ++j) { for (size_t j = 0; j < channels[i].size(); ++j) {
auto &ch = _channels[i][j]; auto &ch = channels[i][j];
buffers[i][j].resize(1024); buffers[i][j].resize(1024);
t->SetBranchAddress((ch + "_Nraw").c_str(), &nraws[i][j]); t->SetBranchAddress((ch + "_Nraw").c_str(), &nraws[i][j]);
t->SetBranchAddress((ch + "_raw").c_str(), &buffers[i][j][0]); t->SetBranchAddress((ch + "_raw").c_str(), &buffers[i][j][0]);
...@@ -90,11 +92,11 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi ...@@ -90,11 +92,11 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi
for (int idx : indices) { for (int idx : indices) {
if (idx >= maxn) { continue; } if (idx >= maxn) { continue; }
t->GetEntry(idx); t->GetEntry(idx);
TestEvent event{"event_" + std::to_string(idx), std::vector<std::vector<TestChannel>>(_channels.size())}; TestEvent event{"event_" + std::to_string(idx), std::vector<std::vector<TestChannel>>(channels.size())};
for (size_t i = 0; i < _channels.size(); ++i) { for (size_t i = 0; i < channels.size(); ++i) {
event.channels[i].resize(_channels[i].size()); event.channels[i].resize(channels[i].size());
for (size_t j = 0; j < _channels[i].size(); ++j) { for (size_t j = 0; j < channels[i].size(); ++j) {
auto &ch = _channels[i][j]; auto &ch = channels[i][j];
auto &nraw = nraws[i][j]; auto &nraw = nraws[i][j];
auto &buf = buffers[i][j]; auto &buf = buffers[i][j];
event.channels[i][j].name = ch; event.channels[i][j].name = ch;
...@@ -108,14 +110,14 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi ...@@ -108,14 +110,14 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi
void draw_events(const std::string &path = "", const std::string &dir = "./plots", void draw_events(const std::string &path = "", const std::string &dir = "./plots",
const std::vector<int> &indices = {}, int nev = 5, bool quadrant = true, const std::vector<int> &indices = {}, int nev = 5,
size_t res = 3, double thres = 10.0) size_t res = 3, double thres = 10.0)
{ {
gSystem->Exec(("mkdir -p " + dir).c_str()); gSystem->Exec(("mkdir -p " + dir).c_str());
gStyle->SetOptStat(0); gStyle->SetOptStat(0);
wfa::Analyzer ana(res, thres); wfa::Analyzer ana(res, thres);
auto events = get_events(path, indices, nev); auto events = get_events(path, indices, nev, quadrant ? _quads : _sums);
for (auto &event : events) { for (auto &event : events) {
auto c = new TCanvas(event.name.c_str(), event.name.c_str(), 1920, 1080); auto c = new TCanvas(event.name.c_str(), event.name.c_str(), 1920, 1080);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment