update wfanalyzer

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:

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