diff --git a/include/WfAnalyzer.h b/include/WfAnalyzer.h index c2d7636cf8ebce73f14cddf95e9d38d5cbecd217..0dedc7b9766a9e020279c03f4aecd85108146ec0 100644 --- a/include/WfAnalyzer.h +++ b/include/WfAnalyzer.h @@ -97,7 +97,7 @@ public: : _thres(threshold), _res(resolution), _npeds(n_ped_samples), _overflow(max_value) {} ~Analyzer() {} - WFData Analyze(const uint32_t *samples, size_t nsamples) + WFData Analyze(const uint32_t *samples, size_t nsamples) const { WFData data; @@ -170,7 +170,7 @@ public: // 1. check the flat-ness of the first few samples, if good, use it as the pedestal // 2. if not, do the same check for the last few samples // 3. if both are bad, use the peaks excluded spectrum (now using TSpectrum) - Pedestal FindPedestal(const std::vector<double> &buffer, const std::vector<Peak> &/*peaks*/) + Pedestal FindPedestal(const std::vector<double> &buffer, const std::vector<Peak> &/*peaks*/) const { Pedestal ped{0., 0.}; // too few samples, use the minimum value as the pedestal @@ -209,6 +209,16 @@ public: return CalcPedestal(&ybuf[ybuf.size()/5], 3*ybuf.size()/5, 1.0, 3); } + double GetThreshold() const { return _thres; } + size_t GetResolution() const { return _res; } + size_t GetNSamplesPed() const { return _npeds; } + uint32_t GetOverflowValue() const { return _overflow; } + + void SetThreshold(double thres) { _thres = thres; } + void SetResolution(size_t res) { _res = res; } + void SetNSamplesPed(size_t npeds) { _npeds = npeds; } + void SetOverflowValue(uint32_t overflow) { _overflow = overflow; } + private: double _thres; size_t _res, _npeds; diff --git a/include/WfGraph.h b/include/WfGraph.h new file mode 100644 index 0000000000000000000000000000000000000000..294ab9b6ac310f53d10b81cb5ca3f9541469c5f8 --- /dev/null +++ b/include/WfGraph.h @@ -0,0 +1,144 @@ +#pragma once + +#include "WfAnalyzer.h" +#include "TCanvas.h" +#include "TGraph.h" +#include "TMultiGraph.h" +#include "TLegend.h" +#include "TSpectrum.h" +#include <algorithm> +#include <string> +#include <vector> +#include <random> + + +namespace wfa { + +struct LegendEntry +{ + TObject *obj; + std::string label, option; +}; + +struct WFGraph +{ + TMultiGraph *mg; + std::vector<TObject*> objs; + std::vector<LegendEntry> entries; + WFData result; +}; + +WFGraph get_waveform_graph(const Analyzer &ana, const std::vector<uint32_t> &samples, bool show_tspec = false) +{ + WFGraph res; + + res.mg = new TMultiGraph(); + + // raw waveform + auto wf = new TGraph(samples.size()); + for (size_t i = 0; i < samples.size(); ++i) { + wf->SetPoint(i, i, samples[i]); + } + wf->SetLineColor(kRed + 1); + wf->SetLineWidth(3); + wf->SetLineStyle(2); + res.mg->Add(wf, "l"); + res.objs.push_back(dynamic_cast<TObject*>(wf)); + res.entries.emplace_back(LegendEntry{wf, "Raw Samples", "l"}); + + // smoothed waveform + auto wf2 = new TGraph(samples.size()); + auto buf = ana.SmoothSpectrum(&samples[0], samples.size(), ana.GetResolution()); + for (size_t i = 0; i < buf.size(); ++i) { + wf2->SetPoint(i, i, buf[i]); + } + wf2->SetLineColor(kBlue + 1); + wf2->SetLineWidth(3); + 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"}); + + // peak finding + res.result = ana.Analyze(&samples[0], samples.size()); + auto ped = res.result.ped; + auto peaks = res.result.peaks; + auto grm_pos = new TGraph(); + auto grm_neg = new TGraph(); + + for (size_t i = 0; i < peaks.size(); ++i) { + // peak amplitudes + auto peak = peaks[i]; + double range = wf->GetYaxis()-> GetXmax() - wf->GetYaxis()-> GetXmin(); + double height = peak.height + ped.mean + (peak.height > 0 ? 1. : -1.5)*range*0.02; + if (peak.height >= 0) { + grm_pos->SetPoint(grm_pos->GetN(), peak.pos, height); + } else { + grm_neg->SetPoint(grm_neg->GetN(), peak.pos, height); + } + + // peak integrals + auto color = i + 2; + 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(3001); + res.mg->Add(grs, "f"); + res.objs.push_back(dynamic_cast<TObject*>(grs)); + if (i == 0) { + res.entries.emplace_back(LegendEntry{grs, "Peak Integrals", "f"}); + } + } + grm_pos->SetMarkerStyle(23); + grm_pos->SetMarkerSize(1.5); + res.mg->Add(grm_pos, "p"); + res.objs.push_back(dynamic_cast<TObject*>(grm_pos)); + grm_neg->SetMarkerStyle(22); + grm_neg->SetMarkerSize(1.5); + res.mg->Add(grm_neg, "p"); + res.objs.push_back(dynamic_cast<TObject*>(grm_neg)); + res.entries.emplace_back(LegendEntry{grm_pos, "Peaks", "p"}); + + // pedestal line + auto grp = new TGraphErrors(samples.size()); + for (size_t i = 0; i < samples.size(); ++i) { + grp->SetPoint(i, i, ped.mean); + grp->SetPointError(i, 0, ped.err); + } + grp->SetFillColor(kBlack); + grp->SetFillStyle(3001); + grp->SetLineStyle(0); + grp->SetLineWidth(2); + grp->SetLineColor(kBlack); + res.mg->Add(grp, "l3"); + res.objs.push_back(dynamic_cast<TObject*>(grp)); + res.entries.emplace_back(LegendEntry{grp, "Pedestal", "lf"}); + + // TSpectrum background + if (show_tspec) { + std::vector<double> tped(samples.size()); + tped.assign(samples.begin(), samples.end()); + TSpectrum s; + s.Background(&tped[0], samples.size(), samples.size()/4, TSpectrum::kBackDecreasingWindow, + TSpectrum::kBackOrder2, false, TSpectrum::kBackSmoothing3, false); + auto grp2 = new TGraph(samples.size()); + for (size_t i = 0; i < samples.size(); ++i) { + grp2->SetPoint(i, i, tped[i]); + } + grp2->SetLineStyle(2); + grp2->SetLineWidth(2); + grp2->SetLineColor(kBlack); + res.mg->Add(grp2, "l"); + res.objs.push_back(dynamic_cast<TObject*>(grp2)); + res.entries.emplace_back(LegendEntry{grp2, "TSpectrum Background", "l"}); + } + + res.mg->GetXaxis()->SetRangeUser(0, samples.size() - 1); + return res; +} + +} // namespace wfa diff --git a/tools/test_wfanalyzer.cxx b/scripts/root/test_wfanalyzer.cxx similarity index 58% rename from tools/test_wfanalyzer.cxx rename to scripts/root/test_wfanalyzer.cxx index 3045bd289a30965a825fe3653c20fb079da954b8..c99f5198ffef19c3dbe656bed6d150f8ed12f63a 100644 --- a/tools/test_wfanalyzer.cxx +++ b/scripts/root/test_wfanalyzer.cxx @@ -1,4 +1,5 @@ #include "WfAnalyzer.h" +#include "WfGraph.h" #include <algorithm> #include <string> #include <vector> @@ -138,104 +139,17 @@ void test_wfanalyzer(const std::string &path = "", const std::vector<int> &indic int count = 0; for (auto &event : events) { - TMultiGraph *mg = new TMultiGraph(); - mg->SetTitle((event.name + "; Sample Number; ADC Values").c_str()); - auto samples = event.raw; + auto graph = get_waveform_graph(ana, event.raw); + graph.mg->SetTitle((event.name + "; Sample Number; ADC Values").c_str()); count ++; c->cd(count); - // waveform - auto wf = new TGraph(samples.size()); - for (size_t i = 0; i < samples.size(); ++i) { - wf->SetPoint(i, i, samples[i]); - } - wf->SetLineColor(kRed + 1); - wf->SetLineWidth(3); - wf->SetLineStyle(2); - mg->Add(wf, "l"); - - // waveform resolution - auto wf2 = new TGraph(samples.size()); - 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]); - } - wf2->SetLineColor(kBlue + 1); - wf2->SetLineWidth(3); - mg->Add(wf2, "l"); - - // peak finding - auto data = ana.Analyze(&samples[0], samples.size()); - auto ped = data.ped; - for (size_t i = 0; i < data.peaks.size(); ++i) { - auto color = i + 2; - auto peak = data.peaks[i]; - auto grm = new TGraph(1); - 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); - mg->Add(grm, "p"); - - 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(3001); - mg->Add(grs, "f"); - // 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"); - } - } + graph.mg->Draw("A"); - // pedestal line - auto grp = new TGraphErrors(samples.size()); - for (size_t i = 0; i < samples.size(); ++i) { - grp->SetPoint(i, i, ped.mean); - grp->SetPointError(i, 0, ped.err); - // grp->SetPoint(i, i, tped[i]); - } - grp->SetFillColor(kBlack); - grp->SetFillStyle(3001); - grp->SetLineStyle(0); - grp->SetLineWidth(2); - grp->SetLineColor(kBlack); - mg->Add(grp, "l3"); - - /* - // TSpectrum background - std::vector<double> tped(samples.size()); - tped.assign(samples.begin(), samples.end()); - TSpectrum s; - s.Background(&tped[0], samples.size(), samples.size()/4, TSpectrum::kBackDecreasingWindow, - TSpectrum::kBackOrder2, false, TSpectrum::kBackSmoothing3, false); - auto grp2 = new TGraph(samples.size()); - for (size_t i = 0; i < samples.size(); ++i) { - grp2->SetPoint(i, i, tped[i]); - } - grp2->SetLineStyle(2); - grp2->SetLineWidth(2); - grp2->SetLineColor(kBlack); - mg->Add(grp2, "l"); - if (count == 1) { legend->AddEntry(grp2, "TSpectrum Background", "l"); } - */ - - if (count == 1) { - legend->AddEntry(wf, "Raw Samples", "l"); - legend->AddEntry(wf2, Form("Smoothed Samples (res = %zu)", res), "l"); - legend->AddEntry(grp, "Pedestal", "lf"); + if (count > 1) { continue; } + + for (auto &entry : graph.entries) { + legend->AddEntry(entry.obj, entry.label.c_str(), entry.option.c_str()); } - mg->GetXaxis()->SetRangeUser(0, samples.size() - 1); - mg->Draw("A"); } c->cd(events.size() + 1); legend->Draw();