diff --git a/include/WfAnalyzer.h b/include/WfAnalyzer.h index fa3b0bf19a8c715d0ab07eddf499e56f638cf5f4..7cfe326defd08849fef50a46acfe57700efc23e5 100644 --- a/include/WfAnalyzer.h +++ b/include/WfAnalyzer.h @@ -111,13 +111,15 @@ public: // get final results for (auto &peak : candidates) { - // real sample height - 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*data.ped.err)) { + // pedestal subtraction + double peak_height = buffer[peak.pos] - data.ped.mean; + // wrong baselin in the rough candidtes finding, below threshold, or not statistically significant + if ((peak_height * peak.height < 0.) || + (std::abs(peak_height) < _thres) || + (std::abs(peak_height) < 3.0*data.ped.err)) { continue; } + peak.height = peak_height; // integrate it over the peak range peak.integral = buffer[peak.pos] - data.ped.mean; @@ -125,28 +127,28 @@ public: for (int i = peak.pos - 1; i >= static_cast<int>(peak.left); --i) { double val = buffer[i] - data.ped.mean; // 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.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 || val * sample_height < 0.) { + if (std::abs(val) < data.ped.err || val * peak.height < 0.) { peak.right = i; break; } peak.integral += val; } // 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) { 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 (sample_pos > i) { update_peak(peak, samples[sample_pos - i] - data.ped.mean, sample_pos - i); @@ -189,7 +191,7 @@ public: // progressively find a good baseline ped.mean = max_mean; 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); if(err < good_err && mean < max_mean) { find_baseline = true; @@ -272,7 +274,7 @@ public: std::vector<Peak> candidates; if (buffer.size() < 3) { return candidates; } - candidates.reserve(buffer.size()/2); + candidates.reserve(buffer.size()/3); // get trend auto trend = [] (double v1, double v2, double thr = 0.1) { return std::abs(v1 - v2) < thr ? 0 : (v1 > v2 ? 1 : -1); @@ -291,8 +293,7 @@ public: right ++; } - double base = buffer[i - left] * right / static_cast<double>(left + right) - + buffer[i + right] * left / static_cast<double>(left + right); + double base = (buffer[i - left] * right + buffer[i + right] * left) / static_cast<double>(left + right); double height = std::abs(buffer[i] - base); if (height > height_thres) { candidates.push_back(Peak{i, i - left, i + right, buffer[i] - base, 0}); diff --git a/include/WfGraph.h b/include/WfGraph.h index e679846b776a490a6b9a210a6d90ca0a1eff25be..d9d807ef226fcb3dd2e8cbf452f25a057c571064 100644 --- a/include/WfGraph.h +++ b/include/WfGraph.h @@ -28,7 +28,8 @@ struct WFGraph 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; @@ -47,7 +48,7 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam wf->SetPoint(i, i, samples[i]); } wf->SetLineColor(kRed + 1); - wf->SetLineWidth(3); + wf->SetLineWidth(1); wf->SetLineStyle(2); res.mg->Add(wf, "l"); 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 wf2->SetPoint(i, i, buf[i]); } wf2->SetLineColor(kBlue + 1); - wf2->SetLineWidth(3); + wf2->SetLineWidth(2); res.mg->Add(wf2, "l"); res.objs.push_back(dynamic_cast<TObject*>(wf2)); 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 grs->SetPoint(i, i + peak.left, val); grs->SetPoint(nint + i, peak.right - i, ped.mean); } - grs->SetFillColor(color); - grs->SetFillStyle(3001); + grs->SetFillColorAlpha(color, shade_alpha); + grs->SetFillStyle(3000); res.mg->Add(grs, "f"); res.objs.push_back(dynamic_cast<TObject*>(grs)); if (i == 0) { @@ -101,13 +102,13 @@ WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &sam } } grm_pos->SetMarkerStyle(23); - grm_pos->SetMarkerSize(1.5); + grm_pos->SetMarkerSize(1.0); if (grm_pos->GetN()) { res.mg->Add(grm_pos, "p"); } res.objs.push_back(dynamic_cast<TObject*>(grm_pos)); res.entries.emplace_back(LegendEntry{grm_pos, "Peaks", "p"}); grm_neg->SetMarkerStyle(22); - grm_neg->SetMarkerSize(1.5); + grm_neg->SetMarkerSize(1.0); if (grm_neg->GetN()) { res.mg->Add(grm_neg, "p"); } 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 grp->SetPoint(i, i, ped.mean); grp->SetPointError(i, 0, ped.err); } - grp->SetFillColor(kBlack); - grp->SetFillStyle(3001); - grp->SetLineStyle(0); + grp->SetFillColorAlpha(kBlack, shade_alpha); + grp->SetFillStyle(3000); grp->SetLineWidth(2); grp->SetLineColor(kBlack); res.mg->Add(grp, "l3"); diff --git a/scripts/root/draw_events.cxx b/scripts/root/draw_events.cxx index 9be0bd942115e1ef88883927523a54a7a6a200a9..b170f03500a86fd541452aa2331fe8ab06151577 100644 --- a/scripts/root/draw_events.cxx +++ b/scripts/root/draw_events.cxx @@ -8,7 +8,7 @@ // 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"}, {"Cer12_1", "Cer12_2", "Cer12_3", "Cer12_4"}, {"Cer13_1", "Cer13_2", "Cer13_3", "Cer13_4"}, @@ -28,14 +28,15 @@ static const std::vector<std::vector<std::string>> _channels = { {"Cer42_1", "Cer42_2", "Cer42_3", "Cer42_4"}, {"Cer43_1", "Cer43_2", "Cer43_3", "Cer43_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 { std::string name; std::vector<uint32_t> raw; @@ -47,7 +48,8 @@ struct TestEvent { }; // 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; // root tree @@ -56,13 +58,13 @@ std::vector<TestEvent> get_events(const std::string &path, std::vector<int> indi int maxn = t->GetEntries(); // fetch tree data - std::vector<std::vector<int>> nraws(_channels.size()); - std::vector<std::vector<std::vector<int>>> buffers(_channels.size()); - for (size_t i = 0; i < _channels.size(); ++i) { - nraws[i].resize(_channels[i].size()); - buffers[i].resize(_channels[i].size()); - for (size_t j = 0; j < _channels[i].size(); ++j) { - auto &ch = _channels[i][j]; + std::vector<std::vector<int>> nraws(channels.size()); + std::vector<std::vector<std::vector<int>>> buffers(channels.size()); + for (size_t i = 0; i < channels.size(); ++i) { + nraws[i].resize(channels[i].size()); + buffers[i].resize(channels[i].size()); + for (size_t j = 0; j < channels[i].size(); ++j) { + auto &ch = channels[i][j]; buffers[i][j].resize(1024); t->SetBranchAddress((ch + "_Nraw").c_str(), &nraws[i][j]); 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 for (int idx : indices) { if (idx >= maxn) { continue; } t->GetEntry(idx); - TestEvent event{"event_" + std::to_string(idx), std::vector<std::vector<TestChannel>>(_channels.size())}; - for (size_t i = 0; i < _channels.size(); ++i) { - event.channels[i].resize(_channels[i].size()); - for (size_t j = 0; j < _channels[i].size(); ++j) { - auto &ch = _channels[i][j]; + TestEvent event{"event_" + std::to_string(idx), std::vector<std::vector<TestChannel>>(channels.size())}; + for (size_t i = 0; i < channels.size(); ++i) { + event.channels[i].resize(channels[i].size()); + for (size_t j = 0; j < channels[i].size(); ++j) { + auto &ch = channels[i][j]; auto &nraw = nraws[i][j]; auto &buf = buffers[i][j]; event.channels[i][j].name = ch; @@ -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", - 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) { gSystem->Exec(("mkdir -p " + dir).c_str()); gStyle->SetOptStat(0); 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) { auto c = new TCanvas(event.name.c_str(), event.name.c_str(), 1920, 1080);