diff --git a/config/#esb_module.conf# b/config/#esb_module.conf#
deleted file mode 100644
index 8c3cbeaef893cf04aa1a8e481615826776ad0180..0000000000000000000000000000000000000000
--- a/config/#esb_module.conf#
+++ /dev/null
@@ -1,208 +0,0 @@
-# configuration of DAQ modules
-# Module { crate, slot, mod_type, Channels { id, name, ch_type } }
-
-Module
-{
- crate = 1
- slot = 3
- type = FADC250
- Channels
- {
- 0, Cer14_5, Signal
- 1, Cer13_4, Signal
- 2, Cer13_3, Signal
- 3, Cer13_2, Signal
- 4, Cer13_1, Signal
- 5, Cer13_5, Signal
- 6, Cer12_4, Signal
- 7, Cer12_3, Signal
- 8, Cer12_2, Signal
- 9, Cer12_1, Signal
- 10, Cer12_5, Signal
- 11, Cer11_4, Signal
- 12, Cer11_3, Signal
- 13, Cer11_2, Signal
- 14, Cer11_1, Signal
- 15, Cer11_5, Signal
- }
-}
-
-Module
-{
- crate = 1
- slot = 4
- type = FADC250
- Channels
- {
- 0, Cer33_1, Signal
- 1, Cer33_5, Signal
- 2, Cer23_4, Signal
- 3, Cer23_3, Signal
- 4, Cer23_2, Signal
- 5, Cer23_1, Signal
- 6, Cer23_5, Signal
- 7, Cer24_4, Signal
- 8, Cer24_3, Signal
- 9, Cer24_2, Signal
- 10, Cer24_1, Signal
- 11, Cer24_5, Signal
- 12, Cer14_4, Signal
- 13, Cer14_3, Signal
- 14, Cer14_2, Signal
- 15, Cer14_1, Signal
- }
-}
-
-Module
-{
- crate = 1
- slot = 5
- type = FADC250
- Channels
- {
- 0, Cer43_2, Signal
- 1, Cer43_1, Signal
- 2, Cer43_5, Signal
- 3, Cer44_4, Signal
- 4, Cer44_3, Signal
- 5, Cer44_2, Signal
- 6, Cer44_1, Signal
- 7, Cer44_5, Signal
- 8, Cer34_4, Signal
- 9, Cer34_3, Signal
- 10, Cer34_2, Signal
- 11, Cer34_1, Signal
- 12, Cer34_5, Signal
- 13, Cer33_4, Signal
- 14, Cer33_3, Signal
- 15, Cer33_2, Signal
- }
-}
-
-Module
-{
- crate = 1
- slot = 6
- type = FADC250
- Channels
- {
- 0, Cer31_3, Signal
- 1, Cer31_2, Signal
- 2, Cer31_1, Signal
- 3, Cer31_5, Signal
- 4, Cer41_4, Signal
- 5, Cer41_3, Signal
- 6, Cer41_2, Signal
- 7, Cer41_1, Signal
- 8, Cer41_5, Signal
- 9, Cer42_4, Signal
- 10, Cer42_3, Signal
- 11, Cer42_2, Signal
- 12, Cer42_1, Signal
- 13, Cer42_5, Signal
- 14, Cer43_4, Signal
- 15, Cer43_3, Signal
- }
-}
-
-Module
-{
- crate = 1
- slot = 7
- type = FADC250
- Channels
- {
- 0, C7_4, Trigger
- 1, C7_3, Trigger
- 2, C7_2, Trigger
- 3, C7_1, Trigger
- 4, C6_4, Trigger
- 5, C6_3, Trigger
- 6, C6_2, Trigger
- 7, C6_1, Trigger
- 8, C5_4, Trigger
- 9, C5_3, Trigger
- 10, C5_2, Trigger
- 11, C5_1, Trigger
- 12, C4, Trigger
- 13, C3, Trigger
- 14, C2, Trigger
- 15, C1, Trigger
- }
-}
-
-Module
-{
- crate = 1
- slot = 8
- type = FADC250
- Channels
- {
- 0, Cer22_1, Signal
- 1, Cer22_5, Signal
- 2, Cer32_4, Signal
- 3, Cer32_3, Signal
- 4, Cer32_2, Signal
- 5, Cer32_1, Signal
- 6, Cer32_5, Signal
- 7, Cer31_4, Signal
- 8, C9_4, Trigger
- 9, C9_3, Trigger
- 10, C9_2, Trigger
- 11, C9_1, Trigger
- 12, C8_4, Trigger
- 13, C8_3, Trigger
- 14, C8_2, Trigger
- 15, C8_1, Trigger
- }
-}
-
-Module
-{
- crate = 1
- slot = 9
- type = FADC250
- Channels
- {
- 2, Cer21_4, Signal
- 3, Cer21_3, Signal
- 4, Cer21_2, Signal
- 5, Cer21_1, Signal
- 6, Cer21_5, Signal
- 7, Cer22_4, Signal
- 8, Cer22_3, Signal
- 9, Cer22_2, Signal
- 10, S2_11, Trigger
- 11, S2_10, Trigger
- 12, S2_9, Trigger
- 13, S2_8, Trigger
- 14, S2_7, Trigger
- 15, S2_6, Trigger
- }
-}
-
-Module
-{
- crate = 1
- slot = 10
- type = FADC250
- Channels
- {
- 0, S2_5, Trigger
- 1, S2_4, Trigger
- 2, S2_3, Trigger
- 3, S2_2, Trigger
- 4, S2_1, Trigger
- 5, S1_11, Trigger
- 6, S1_10, Trigger
- 7, S1_9, Trigger
- 8, S1_8, Trigger
- 9, S1_7, Trigger
- 10, S1_6, Trigger
- 11, S1_5, Trigger
- 12, S1_4, Trigger
- 13, S1_3, Trigger
- 14, S1_2, Trigger
- 15, S1_1, Trigger
- }
-}
diff --git a/config/esb_module.conf.bkp b/config/esb_module_old.conf
similarity index 100%
rename from config/esb_module.conf.bkp
rename to config/esb_module_old.conf
diff --git a/config/json/esb_module.json b/config/json/esb_module.json
new file mode 100644
index 0000000000000000000000000000000000000000..ef72db1c557ce50552854961f498c9014480ce5b
--- /dev/null
+++ b/config/json/esb_module.json
@@ -0,0 +1,688 @@
+{
+ "module_1_10": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "S2_5",
+ "type": "Trigger"
+ },
+ {
+ "channel": 1,
+ "name": "S2_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 2,
+ "name": "S2_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 3,
+ "name": "S2_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 4,
+ "name": "S2_1",
+ "type": "Trigger"
+ },
+ {
+ "channel": 5,
+ "name": "S1_11",
+ "type": "Trigger"
+ },
+ {
+ "channel": 6,
+ "name": "S1_10",
+ "type": "Trigger"
+ },
+ {
+ "channel": 7,
+ "name": "S1_9",
+ "type": "Trigger"
+ },
+ {
+ "channel": 8,
+ "name": "S1_8",
+ "type": "Trigger"
+ },
+ {
+ "channel": 9,
+ "name": "S1_7",
+ "type": "Trigger"
+ },
+ {
+ "channel": 10,
+ "name": "S1_6",
+ "type": "Trigger"
+ },
+ {
+ "channel": 11,
+ "name": "S1_5",
+ "type": "Trigger"
+ },
+ {
+ "channel": 12,
+ "name": "S1_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 13,
+ "name": "S1_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 14,
+ "name": "S1_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 15,
+ "name": "S1_1",
+ "type": "Trigger"
+ }
+ ],
+ "crate": 1,
+ "slot": 10,
+ "type": "FADC250"
+ },
+ "module_1_3": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "Cer14_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 1,
+ "name": "Cer13_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer13_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer13_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer13_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer13_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer12_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer12_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "Cer12_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 9,
+ "name": "Cer12_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 10,
+ "name": "Cer12_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 11,
+ "name": "Cer11_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 12,
+ "name": "Cer11_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 13,
+ "name": "Cer11_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 14,
+ "name": "Cer11_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 15,
+ "name": "Cer11_5",
+ "type": "Signal"
+ }
+ ],
+ "crate": 1,
+ "slot": 3,
+ "type": "FADC250"
+ },
+ "module_1_4": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "Cer33_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 1,
+ "name": "Cer33_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer23_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer23_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer23_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer23_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer23_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer24_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "Cer24_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 9,
+ "name": "Cer24_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 10,
+ "name": "Cer24_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 11,
+ "name": "Cer24_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 12,
+ "name": "Cer14_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 13,
+ "name": "Cer14_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 14,
+ "name": "Cer14_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 15,
+ "name": "Cer14_1",
+ "type": "Signal"
+ }
+ ],
+ "crate": 1,
+ "slot": 4,
+ "type": "FADC250"
+ },
+ "module_1_5": {
+ "channels": [
+ {
+ "channel": 1,
+ "name": "Cer43_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer43_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer44_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer44_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer44_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer44_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer44_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "Cer34_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 9,
+ "name": "Cer34_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 10,
+ "name": "Cer34_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 11,
+ "name": "Cer34_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 12,
+ "name": "Cer34_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 13,
+ "name": "Cer33_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 14,
+ "name": "Cer33_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 15,
+ "name": "Cer33_2",
+ "type": "Signal"
+ }
+ ],
+ "crate": 1,
+ "slot": 5,
+ "type": "FADC250"
+ },
+ "module_1_6": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "Cer31_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 1,
+ "name": "Cer31_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer31_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer31_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer41_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer41_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer41_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer41_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "Cer41_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 9,
+ "name": "Cer42_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 10,
+ "name": "Cer42_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 11,
+ "name": "Cer42_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 12,
+ "name": "Cer42_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 13,
+ "name": "Cer42_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 14,
+ "name": "Cer43_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 15,
+ "name": "Cer43_2",
+ "type": "Signal"
+ }
+ ],
+ "crate": 1,
+ "slot": 6,
+ "type": "FADC250"
+ },
+ "module_1_7": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "C7_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 1,
+ "name": "C7_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 2,
+ "name": "C7_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 3,
+ "name": "C7_1",
+ "type": "Trigger"
+ },
+ {
+ "channel": 4,
+ "name": "C6_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 5,
+ "name": "C6_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 6,
+ "name": "C6_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 7,
+ "name": "C6_1",
+ "type": "Trigger"
+ },
+ {
+ "channel": 8,
+ "name": "C5_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 9,
+ "name": "C5_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 10,
+ "name": "C5_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 11,
+ "name": "C5_1",
+ "type": "Trigger"
+ },
+ {
+ "channel": 12,
+ "name": "C4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 13,
+ "name": "C3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 14,
+ "name": "C2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 15,
+ "name": "C1",
+ "type": "Trigger"
+ }
+ ],
+ "crate": 1,
+ "slot": 7,
+ "type": "FADC250"
+ },
+ "module_1_8": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "Cer22_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 1,
+ "name": "Cer22_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer32_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer32_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer32_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer32_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer32_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer31_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "C9_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 9,
+ "name": "C9_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 10,
+ "name": "C9_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 11,
+ "name": "C9_1",
+ "type": "Trigger"
+ },
+ {
+ "channel": 12,
+ "name": "C8_4",
+ "type": "Trigger"
+ },
+ {
+ "channel": 13,
+ "name": "C8_3",
+ "type": "Trigger"
+ },
+ {
+ "channel": 14,
+ "name": "C8_2",
+ "type": "Trigger"
+ },
+ {
+ "channel": 15,
+ "name": "C8_1",
+ "type": "Trigger"
+ }
+ ],
+ "crate": 1,
+ "slot": 8,
+ "type": "FADC250"
+ },
+ "module_1_9": {
+ "channels": [
+ {
+ "channel": 0,
+ "name": "Cer43_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 2,
+ "name": "Cer21_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 3,
+ "name": "Cer21_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 4,
+ "name": "Cer21_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 5,
+ "name": "Cer21_1",
+ "type": "Signal"
+ },
+ {
+ "channel": 6,
+ "name": "Cer21_5",
+ "type": "Signal"
+ },
+ {
+ "channel": 7,
+ "name": "Cer22_4",
+ "type": "Signal"
+ },
+ {
+ "channel": 8,
+ "name": "Cer22_3",
+ "type": "Signal"
+ },
+ {
+ "channel": 9,
+ "name": "Cer22_2",
+ "type": "Signal"
+ },
+ {
+ "channel": 10,
+ "name": "S2_11",
+ "type": "Trigger"
+ },
+ {
+ "channel": 11,
+ "name": "S2_10",
+ "type": "Trigger"
+ },
+ {
+ "channel": 12,
+ "name": "S2_9",
+ "type": "Trigger"
+ },
+ {
+ "channel": 13,
+ "name": "S2_8",
+ "type": "Trigger"
+ },
+ {
+ "channel": 14,
+ "name": "S2_7",
+ "type": "Trigger"
+ },
+ {
+ "channel": 15,
+ "name": "S2_6",
+ "type": "Trigger"
+ }
+ ],
+ "crate": 1,
+ "slot": 9,
+ "type": "FADC250"
+ }
+}
diff --git a/include/WfAnalyzer.h b/include/WfAnalyzer.h
new file mode 100644
index 0000000000000000000000000000000000000000..c6f5b7f4a0c902138288d3523dbce6b5f04c215b
--- /dev/null
+++ b/include/WfAnalyzer.h
@@ -0,0 +1,242 @@
+#pragma once
+
+// A waveform analyzer class for the Fadc250 raw waveform data
+// It finds the local maxima (peaks) over the spectrum, and detemines the peak properties
+// It is intended to use simple algorithms since we do not expect heavily convoluted spectra
+//
+// History:
+// v1.0.0: a working version with basic functionalities, Chao Peng, 2020/08/04
+
+#include <vector>
+#include <algorithm>
+
+// calculate the mean and standard deviation for an array
+namespace wfa {
+
+// data structures
+struct Pedestal {
+ double mean, err;
+};
+
+struct Peak {
+ size_t pos, left, right;
+ double height, integral;
+
+ bool Inside(size_t i) {
+ return (i >= pos - left) && (i <= pos + right);
+ }
+};
+
+// 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)
+{
+ if (npts == 0) { return; }
+
+ mean = 0.;
+ err = 0.;
+
+ for (size_t i = 0; i < npts; ++i) {
+ mean += source[i];
+ }
+ mean /= static_cast<double>(npts);
+
+ for (size_t i = 0; i < npts; ++i) {
+ err += (source[i] - mean ) * (source[i] - mean);
+ }
+ err = std::sqrt(err/static_cast<double>(npts));
+}
+
+template<typename T>
+inline void _linear_regr(double &p0, double &p1, double &err, T *x, T *y, size_t npts)
+{
+ // no data points
+ if (npts == 0) {
+ return;
+ // too few data points
+ } else if (npts < 3) {
+ p1 = 0.;
+ _calc_mean_err(p0, err, y, npts);
+ return;
+ }
+
+ double x_mean, x_err;
+ _calc_mean_err(x_mean, x_err, x, npts);
+ double y_mean, y_err;
+ _calc_mean_err(y_mean, y_err, y, npts);
+
+ double xysum = 0.;
+ for (size_t i = 0; i < npts; ++i) {
+ xysum += (x[i] - x_mean) * (y[i] - y_mean);
+ }
+
+ double r = xysum/(x_err*y_err*npts);
+
+ p1 = r*y_err/x_err;
+ p0 = y_mean - p1*x_mean;
+
+ err = 0.;
+ for (size_t i = 0; i < npts; ++i) {
+ double diff = (p0 + p1*x[i] - y[i]);
+ err += diff*diff;
+ }
+ err = std::sqrt(err/static_cast<double>(npts));
+}
+
+// analyzer class
+class Analyzer
+{
+public:
+ Analyzer(size_t resolution = 2, double threshold = 10.0)
+ : thres(threshold), res(resolution) {}
+ ~Analyzer() {}
+
+ void Analyze(uint32_t *samples, size_t nsamples)
+ {
+ peaks.clear();
+ if (!nsamples) { return; }
+
+ auto buffer = ResSpectrum(samples, nsamples, res);
+
+ // search local maxima
+ std::vector<Peak> candidates;
+ candidates.reserve(nsamples/2);
+ // get trend
+ auto trend = [] (double v1, double v2, double thr = 0.1) {
+ return std::abs(v1 - v2) < thr ? 0 : (v1 > v2 ? 1 : -1);
+ };
+ for (size_t i = 1; i < nsamples - 1; ++i) {
+ 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) {
+ size_t left = 1, right = 1;
+ // search the peak range
+ while ((i > left + 1) && (trend(buffer[i - left], buffer[i - left - 1]) == tr1)) {
+ left ++;
+ }
+ while ((i + right < nsamples - 1) && (trend(buffer[i + right], buffer[i + right + 1])*tr1 >= 0)) {
+ right ++;
+ }
+
+ double base = buffer[i - left] * right / static_cast<double>(left + right)
+ + buffer[i + right] * left / static_cast<double>(left + right);
+ double height = std::abs(buffer[i] - base);
+ if (height > thres) {
+ candidates.push_back(Peak{i, left, right, buffer[i] - base, 0});
+ }
+ }
+ }
+
+ auto in_peaks = [] (const std::vector<Peak> &peaks, size_t pos) {
+ for (auto &peak : peaks) {
+ size_t width = std::min(peak.left, peak.right);
+ if (std::abs(peak.pos - static_cast<double>(pos)) <= width) { return true; }
+ }
+ return false;
+ };
+
+ int nped = 0;
+ auto ybuf = buffer;
+ for (size_t i = res + 1; i < nsamples - res - 1; ++i) {
+ if (!in_peaks(candidates, i)) {
+ ybuf[nped] = buffer[i];
+ nped ++;
+ }
+ }
+
+ // calculate pedestal
+ if (nped <= nsamples/4) {
+ ped = CalcPedestal(&buffer[0], nsamples, 1.5, 3);
+ } else {
+ ped = CalcPedestal(&ybuf[0], nped, 1.2, 3);
+ }
+
+ for (auto &peak : candidates) {
+ // real sample height
+ double sample_height = samples[peak.pos] - ped.mean;
+ if ((sample_height * peak.height < 0.) ||
+ (std::abs(sample_height) < thres) ||
+ (std::abs(sample_height) < 3.0*ped.err)) {
+ continue;
+ }
+
+ // integrate it over the peak range
+ peak.height = sample_height;
+ peak.integral = buffer[peak.pos] - 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; }
+ 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; }
+ peak.integral += val;
+ }
+ peaks.emplace_back(peak);
+ }
+
+ /*
+ std::sort(peaks.begin(), peaks.end(),
+ [] (const Peak &p1, const Peak &p2) { return p1.height > p2.height; });
+ */
+
+ }
+
+ 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:
+ // static methods
+ template<typename T>
+ static std::vector<double> ResSpectrum(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) {
+ if (j >= i || j + i >= nsamples) { continue; }
+ double weight = j/static_cast<double>(res + 1);
+ val += weight*(samples[i - j] + samples[i + j]);
+ weights += 2.*weight;
+ }
+ buffer[i] = val/weights;
+ }
+ return buffer;
+ }
+
+ template<typename T>
+ static Pedestal CalcPedestal(T *ybuf, size_t npts, double thres = 1.0, int max_iters = 3)
+ {
+ Pedestal res;
+ _calc_mean_err(res.mean, res.err, ybuf, npts);
+ // interatively fit
+ while (max_iters-- > 0) {
+ int count = 0;
+ for (size_t i = 0; i < npts; ++i) {
+ if (std::abs(ybuf[i] - res.mean) < res.err*thres) {
+ ybuf[count] = ybuf[i];
+ count++;
+ }
+ }
+ if (count == npts) { break; }
+ npts = count;
+ _calc_mean_err(res.mean, res.err, ybuf, npts);
+ }
+
+ return res;
+ }
+
+};
+
+}; // namespace wfa
+
diff --git a/src/esb_analyze_online.cpp b/src/esb_analyze_online.cpp
deleted file mode 100644
index f7a896264a78ac536f2c11e468ec183ce417df22..0000000000000000000000000000000000000000
--- a/src/esb_analyze_online.cpp
+++ /dev/null
@@ -1,433 +0,0 @@
-#include <iostream>
-#include <iomanip>
-#include <fstream>
-#include <csignal>
-#include "utils.h"
-#include "ConfigParser.h"
-#include "ConfigObject.h"
-#include "ConfigOption.h"
-#include "TSpectrum.h"
-#include "TTree.h"
-#include "TFile.h"
-#include "TH1.h"
-#include "simpleLib.h"
-#include "Fadc250Decoder.h"
-#include "ETChannel.h"
-
-
-#define FADC_BANK 3
-#define PROGRESS_COUNT 100
-
-
-volatile sig_atomic_t terminate = 0;
-
-void handle_sig(int signum)
-{
- terminate = 1;
-}
-
-
-
-void monitor(const std::string &addr, int port, const std::string &etfile,
- const std::string &opath, const std::string &modules_conf, int nev);
-
-
-uint32_t swap_endian32(uint32_t num)
-{
- uint32_t b0 = (num & 0x000000ff) << 24u;
- uint32_t b1 = (num & 0x0000ff00) << 8u;
- uint32_t b2 = (num & 0x00ff0000) >> 8u;
- uint32_t b3 = (num & 0xff000000) >> 24u;
- return b0 | b1 | b2 | b3;
-}
-
-
-uint32_t parseBlock(const uint32_t *buf)
-{
- auto header = (ETChannel::CodaEvHeader*) buf;
- const uint32_t buf_size = header->length + 1;
- std::cout << "ROC header: " << std::dec << buf_size << "\n"
- << (int) header->etype << ", " << (int) header->dtype << ", " << (int) header->num
- << std::endl;
- std::cout << std::hex;
- for (uint32_t i = 0; i < buf_size; ++i) {
- std::cout << "0x" << std::setw(8) << std::setfill('0') << buf[i] << "\n";
- }
- return buf_size;
-}
-
-// parse an evio event
-bool parseEvent(const uint32_t *buf, bool verbose = false)
-{
- auto header = (ETChannel::CodaEvHeader*)buf;
- const uint32_t buf_size = header->length + 1;
- if (verbose) {
- std::cout << "Event header: " << std::dec << buf_size << "\n"
- << (int) header->etype << ", " << (int) header->dtype << ", " << (int) header->num
- << std::endl;
- std::cout << std::hex;
- for (uint32_t i = 0; i < 2; ++i) {
- std::cout << "0x" << std::setw(8) << std::setfill('0') << buf[i] << "\n";
- }
- }
-
- switch(header->etype) {
- case ETChannel::CODA_PHY1:
- case ETChannel::CODA_PHY2:
- break;
- case ETChannel::CODA_PRST:
- case ETChannel::CODA_GO:
- case ETChannel::CODA_END:
- default:
- return false;
- }
-
- simpleScan((volatile uint32_t*)buf, buf_size);
- return true;
-
- /*
- // parse ROC data
- uint32_t index = 2;
- while(index < buf_size) {
- // skip header size and data size 2 + (length - 1)
- index += parseBlock(&buf[index]);
- }
-
- // return parsed event type
- return buf_size;
- */
-}
-
-
-int main(int argc, char* argv[])
-{
- // setup input arguments
- ConfigOption conf_opt;
- conf_opt.AddOpts(ConfigOption::help_message, 'h', "help");
- conf_opt.AddOpt(ConfigOption::arg_require, 'n');
- conf_opt.AddOpt(ConfigOption::arg_require, 'p');
- conf_opt.AddOpt(ConfigOption::arg_require, 'f');
- conf_opt.AddLongOpt(ConfigOption::arg_require, "config-module", 'm');
- conf_opt.AddLongOpt(ConfigOption::arg_require, "host", 'o');
-
- conf_opt.SetDesc("usage: %0 <data_file> <out_file>");
- conf_opt.SetDesc('n', "number of events to process (< 0 means all).");
- conf_opt.SetDesc('m', "configuration file for modules to be read-in, default \"config/esb_module.conf\".");
- conf_opt.SetDesc('p', "port to et system, default 11111.");
- conf_opt.SetDesc('f', "path to et file, default \"/tmp/et_sys\"");
- conf_opt.SetDesc('o', "host address of et syste, default localhost");
-
- if (!conf_opt.ParseArgs(argc, argv) || conf_opt.NbofArgs() != 1) {
- std::cout << conf_opt.GetInstruction() << std::endl;
- return -1;
- }
-
- int nev = -1;
- std::string mconf = "config/esb_module.conf";
- std::string host = "localhost";
- int port = 11111;
- std::string etfile = "/tmp/et_sys";
-
- for (auto &opt : conf_opt.GetOptions()) {
- switch (opt.mark) {
- case 'n':
- nev = opt.var.Int();
- break;
- case 'm':
- mconf = opt.var.String();
- break;
- case 'o':
- host = opt.var.String();
- break;
- case 'p':
- port = opt.var.Int();
- break;
- case 'f':
- etfile = opt.var.String();
- break;
- default :
- std::cout << conf_opt.GetInstruction() << std::endl;
- return -1;
- }
- }
-
- monitor(host, port, etfile, conf_opt.GetArgument(0), mconf, nev);
-}
-
-
-#define MAX_NPEAKS 20
-#define MAX_RAW 300
-struct BranchData
-{
- int npul, nraw;
- float integral[MAX_NPEAKS], peak[MAX_NPEAKS], time[MAX_NPEAKS];
- int raw[MAX_RAW];
- float ped_mean, ped_err;
-};
-
-#define BUF_SIZE 1000
-static double buffer[BUF_SIZE], wfbuf[BUF_SIZE], bkbuf[BUF_SIZE];
-void refine_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
-{
- int count = 0;
- float mean = 0.;
- for (auto &val : raw) {
- if (std::abs(val - ped_mean) < 2.0 * ped_err) {
- buffer[count] = val;
- mean += val;
- count ++;
- }
- }
- if (count == 0) {
- return;
- }
-
- mean /= count;
- float change = std::abs(ped_mean - mean);
-
- // no outliers
- if (change < 0.05 * ped_err) {
- return;
- }
-
- // recursively refine
- float err = 0.;
- for (int i = 0; i < count; ++i) {
- err += (buffer[i] - mean)*(buffer[i] - mean);
- }
- ped_mean = mean;
- ped_err = std::sqrt(err/count);
- refine_pedestal(raw, ped_mean, ped_err);
-}
-
-// a function for a simple constant baseline fit
-void find_pedestal(const std::vector<uint32_t> &raw, float &ped_mean, float &ped_err)
-{
- ped_mean = 0;
- ped_err = 0;
-
- for (auto &val : raw) {
- ped_mean += val;
- }
- ped_mean /= raw.size();
-
- for (auto &val : raw) {
- ped_err += (val - ped_mean)*(val - ped_mean);
- }
- ped_err = std::sqrt(ped_err/raw.size());
-
- refine_pedestal(raw, ped_mean, ped_err);
-}
-
-// analyze the waveform data and fill the result in a branch data
-void waveform_analysis(const std::vector<uint32_t> &raw, BranchData &res)
-{
- // no need to analyze
- if (raw.empty()) {
- return;
- }
-
- // fill in the raw data
- res.nraw = raw.size();
- for (size_t i = 0; i < raw.size(); ++i) {
- res.raw[i] = raw[i];
- }
-
- // get pedestals
- find_pedestal(raw, res.ped_mean, res.ped_err);
-
- // fill in spectrum buffer
- for (size_t i = 0; i < raw.size(); ++i) {
- wfbuf[i] = raw[i] - res.ped_mean;
- }
-
- // find peaks
- TSpectrum s;
- s.SetResolution(0.5);
- int npeaks = s.SearchHighRes(wfbuf, bkbuf, res.nraw, 3.0, 20, false, 5, true, 3);
-
- // fill branch data
- double *pos = s.GetPositionX();
- res.npul = 0;
- for (int i = 0; i < npeaks; ++i) {
- int j = pos[i];
- if (wfbuf[j] < 5.*res.ped_err) { continue; }
- res.time[res.npul] = j;
- res.peak[res.npul] = wfbuf[j];
- /*
- res.integral[res.npul] = wfbuf[j];
- j = pos[i] - 1;
- while ( (j > 0) && (wfbuf[j] - wfbuf[j - 1])*wfbuf[j] > 0. ) {
- res.integral[res.npul] += wfbuf[j--];
- }
- j = pos[i] + 1;
- while ( (j < res.nraw - 1) && (wfbuf[j] - wfbuf[j + 1])*wfbuf[j] > 0. ) {
- res.integral[res.npul] += wfbuf[j++];
- }
- */
- res.npul ++;
- }
-}
-
-
-// fill decoded FADC250 data into the container (branch data)
-void fill_branch(const Fadc250Data &slot_data, const Module &mod, std::unordered_map<std::string, BranchData> &brdata)
-{
- for (auto &event : slot_data.events) {
- for (auto &ch : mod.channels) {
- auto &channel = event.channels[ch.id];
- auto &bd = brdata[ch.name];
- bd.npul = channel.integral.size();
- for (uint32_t i = 0; i < channel.integral.size(); ++i) {
- bd.integral[i] = channel.integral[i];
- bd.time[i] = channel.time[i];
- }
- waveform_analysis(channel.raw, bd);
- }
- }
-}
-
-// fill branch data into the root tree
-void fill_tree(TTree *tree, std::unordered_map<std::string, BranchData> &brdata, bool &init, int mode)
-{
- if ((mode != 1) && (mode != 3)) {
- std::cout << "Warning: unsupported mode " << mode << ", data won't be recorded." << std::endl;
- return;
- }
-
- // initialize
- if (!init) {
- for (auto &it : brdata) {
- auto n = it.first;
- tree->Branch((n + "_Npulse").c_str(), &brdata[n].npul, (n + "_N/I").c_str());
- tree->Branch((n + "_Pint").c_str(), &brdata[n].integral[0], (n + "_Pint[" + n + "_N]/F").c_str());
- tree->Branch((n + "_Ptime").c_str(), &brdata[n].time[0], (n + "_Ptime[" + n + "_N]/F").c_str());
- // raw waveform provides more information
- if (mode == 1) {
- tree->Branch((n + "_Ppeak").c_str(), &brdata[n].peak[0], (n + "_Ppeak[" + n + "_N]/F").c_str());
- tree->Branch((n + "_Nraw").c_str(), &brdata[n].nraw, (n + "_Nraw/I").c_str());
- tree->Branch((n + "_raw").c_str(), &brdata[n].raw[0], (n + "_raw[" + n + "_Nraw]/I").c_str());
- tree->Branch((n + "_ped_mean").c_str(), &brdata[n].ped_mean, (n + "_ped_mean/F").c_str());
- tree->Branch((n + "_ped_err").c_str(), &brdata[n].ped_err, (n + "_ped_err/F").c_str());
- }
- }
- init = true;
- std::cout << "Initialized root file for mode " << mode << " data." << std::endl;
- }
-
- tree->Fill();
-}
-
-
-void processEvent(const uint32_t *buf, int &count, TTree *tree, bool &init_tree, int &data_mode,
- const std::vector<Module> &modules, std::unordered_map<std::string, BranchData> &brdata)
-{
- if (!parseEvent(buf, false)) {
- return;
- }
-
- // get block level
- int blvl = 1;
- simpleGetRocBlockLevel(modules.front().crate, FADC_BANK, &blvl);
-
- for (int ii = 0; ii < blvl; ++ii) {
- // clear data buffer
- for (auto &br : brdata) {
- br.second.npul = 0;
- br.second.nraw = 0;
- }
-
- // parse module data
- for (auto &mod : modules) {
- uint32_t header = 0;
- auto status = simpleGetSlotBlockHeader(mod.crate, FADC_BANK, mod.slot, &header);
- if (status <= 0) {
- std::cout << "Error getting header for crate = "
- << mod.crate << ", slot = " << mod.slot << std::endl;
- continue;
- }
-
- uint32_t *dbuf;
- auto len = simpleGetSlotEventData(mod.crate, FADC_BANK, mod.slot, ii, &dbuf);
- if (len <= 0) {
- std::cout << "No data for crate = " << mod.crate << ", slot = " << mod.slot
- << ", block_level = " << ii << std::endl;
- continue;
- }
- auto slot_data = fadc250_decoder(header, dbuf, len);
-
- // check data mode
- if (data_mode < 0 && slot_data.GetMode() > 0) {
- data_mode = slot_data.GetMode();
- }
- // fill branch data
- fill_branch(slot_data, mod, brdata);
- }
-
- // fill event
- fill_tree(tree, brdata, init_tree, data_mode);
- count ++;
- }
-}
-
-
-void monitor(const std::string &addr, int port, const std::string &etfile,
- const std::string &opath, const std::string &mconf, int nev)
-{
- ETChannel et_chan(100, 100);
- if (et_chan.Connect(addr.c_str(), port, etfile.c_str(), "TCD_MONITOR") != ET_OK) {
- return;
- }
-
- auto modules = read_modules(mconf);
-
- auto *hfile = new TFile(opath.c_str(), "RECREATE", "MAPMT online monitor");
- auto tree = new TTree("EvTree", "Cherenkov Test Events");
- // set up branches and a container to be combined with the branches
- std::unordered_map<std::string, BranchData> brdata;
- for (auto &mod : modules) {
- for (auto &ch : mod.channels) {
- auto n = ch.name;
- if (brdata.find(n) != brdata.end()) {
- std::cout << "WARNING: Duplicated channel names found! -- " << n << std::endl;
- }
- brdata[n] = BranchData();
- }
- }
-
- int count = 0;
- bool init_tree = false;
- int data_mode = -1;
- int status = ETChannel::READ_EMPTY;
- et_chan.Open();
- signal (SIGINT, handle_sig);
-
- while (true) {
- if (terminate > 0 || (nev > 0 && nev < count)) { break; }
- if (count % PROGRESS_COUNT == 0) {
- std::cout << "Read and processed events - " << count << std::endl;
- }
-
- status = et_chan.ReadEvent();
-
- switch (status) {
- case ETChannel::READ_ERROR:
- case ETChannel::READ_EOF:
- terminate = 1;
- break;
- case ETChannel::READ_EMPTY:
- std::cout << "ET station is empty, wait 2 secs..." << std::endl;
- std::this_thread::sleep_for(std::chrono::seconds(2));
- break;
- case ETChannel::READ_OK:
- processEvent(et_chan.GetEvBuffer(), count, tree, init_tree, data_mode, modules, brdata);
- break;
- }
- }
-
- std::cout << "Read and processed events - " << count << std::endl;
-
- hfile->Write();
- hfile->Close();
-}
-
diff --git a/tools/test_wfanalyzer.cxx b/tools/test_wfanalyzer.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..046ee8279a87668728871d017a0499f66b61fb4e
--- /dev/null
+++ b/tools/test_wfanalyzer.cxx
@@ -0,0 +1,95 @@
+#include "WfAnalyzer.h"
+
+
+void test(size_t res = 3, double thres = 10.0)
+{
+ gStyle->SetOptStat(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,
+ 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
+ },
+ {
+ 567, 484, 414, 359, 320, 291, 271, 275, 279, 278, 271, 271, 277,
+ 277, 265, 252, 237, 228, 213, 207, 201, 195, 188, 183, 179, 176,
+ 173, 191, 217, 231, 298, 392, 436, 423, 383, 341, 302, 265, 237,
+ 215, 198, 181, 171, 160, 151, 152, 144, 137, 136, 134, 133, 130,
+ 125, 125, 124, 123, 120, 122, 117, 117, 116, 115, 116, 113
+ },
+ };
+
+ wfa::Analyzer ana(res, thres);
+
+ auto c = new TCanvas("test", "test", 1280, 720);
+ c->DivideSquare(tests.size());
+
+ int count = 0;
+ for (auto &samples : tests) {
+ 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(kBlue + 1);
+ wf->SetLineWidth(1);
+ wf->SetLineStyle(2);
+ wf->Draw("AL");
+
+ // waveform resolution
+ auto wf2 = new TGraph(samples.size());
+ auto buf = wfa::Analyzer::ResSpectrum(&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);
+ wf2->Draw("L same");
+
+ // peak finding
+ ana.Analyze(&samples[0], samples.size());
+ auto ped = ana.GetPedestal();
+ auto gr = new TGraph(samples.size());
+ gr->SetMarkerStyle(23);
+ for (auto &peak : ana.GetPeaks()) {
+ gr->SetPoint(gr->GetN(), peak.pos, peak.height + 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
+ << ", " << peak.height << ", " << peak.integral << std::endl;
+ }
+ gr->Draw("P same");
+
+ // pedestal line
+ /*
+ 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 grp = new TGraph(samples.size());
+ auto gre = new TGraph(2*samples.size());
+ for (size_t i = 0; i < samples.size(); ++i) {
+ double pedy = ped.mean;
+ grp->SetPoint(i, i, pedy);
+ // grp->SetPoint(i, i, tped[i]);
+ gre->SetPoint(i, i, pedy + ped.err);
+ size_t j = samples.size() - i - 1;
+ double pedy2 = ped.mean; // p0 + ped.p1*j;
+ gre->SetPoint(i + samples.size(), j, pedy2 - ped.err);
+ }
+ gre->SetFillStyle(3001);
+ gre->SetFillColor(kRed);
+ grp->SetLineStyle(0);
+ grp->SetLineWidth(2);
+ grp->SetLineColor(kRed);
+ gre->Draw("f same");
+ grp->Draw("L same");
+ }
+}