add LED samples extraction

config/lappd_module.conf

@@ -39,7 +39,7 @@ Module
         2, Calo4, Trigger
         3, Calo5, Trigger
         4, CaloSum, Trigger
-        7, LED, Trigger
+        8, LED, Trigger
         11, Ref, Reference
         12, ScintLR, Trigger
         13, ScintLL, Trigger

config/mapmt_module.conf

@@ -39,7 +39,7 @@ Module
         2, Calo4, Trigger
         3, Calo5, Trigger
         4, CaloSum, Trigger
-        7, LED, Trigger
+        8, LED, Trigger
         11, Ref, Reference
         12, ScintLR, Trigger
         13, ScintLL, Trigger

scripts/batch_run.py

@@ -10,7 +10,8 @@ exec_path = os.path.join(root_path, 'bin/analyze')
 
 # options
 raw_data_fmt = os.path.join(root_path, 'raw_data/fadc_%i.dat')
-output_fmt = os.path.join(root_path, 'processed_data/fadc_%i_proc.root')
+# output_fmt = os.path.join(root_path, 'processed_data/fadc_%i_proc.root')
+output_fmt = os.path.join(root_path, 'raw_data/fadc_%i.root')
 tcut = 2.0
 fire_thres = 30.0
 modules = os.path.join(root_path, 'config/mapmt_module.conf')
@@ -21,15 +22,14 @@ timing = os.path.join(root_path, 'config/mapmt_timing.conf')
 run_options = [
     '--config-module=%s' % modules,
     '--config-timing=%s' % timing,
-    '-c', str(tcut),
-    '-f', str(fire_thres),
-    # '--raw-output'
+    # '-c', str(tcut),
+    # '-f', str(fire_thres),
+    '--raw-output'
 ]
 
 # runs
 # runs = [ 579, 580, 581, 582, 583 ]  # LAPPD
-# runs = [459, 460, 461, 462, 464, 466, 467, 468, 469, 470, 475]  # MAPMT
-runs = [467, 468, 469, 470, 475]
+runs = [459, 460, 461, 462, 464, 466, 467, 468, 469, 470, 475]  # MAPMT
 
 # download setting
 download_host = 'agave.phy.anl.gov'

scripts/draw_samples.py

@@ -9,20 +9,20 @@ import numpy as np
 from matplotlib import pyplot as plt
 from scipy import signal
 
-# from scipy import sparse
-# from scipy.sparse.linalg import spsolve
+from scipy import sparse
+from scipy.sparse.linalg import spsolve
 
 
-# def baseline_als(y, lam=100., p=0.001, niter=10):
-#     L = len(y)
-#     D = sparse.diags([1, -2, 1],[0, -1, -2], shape=(L, L-2))
-#     w = np.ones(L)
-#     for i in range(niter):
-#         W = sparse.spdiags(w, 0, L, L)
-#         Z = W + lam * D.dot(D.transpose())
-#         z = spsolve(Z, w*y)
-#         w = p * (y > z) + (1-p) * (y < z)
-#     return z
+def baseline_als(y, lam=100., p=0.001, niter=10):
+    L = len(y)
+    D = sparse.diags([1, -2, 1],[0, -1, -2], shape=(L, L-2))
+    w = np.ones(L)
+    for i in range(niter):
+        W = sparse.spdiags(w, 0, L, L)
+        Z = W + lam * D.dot(D.transpose())
+        z = spsolve(Z, w*y)
+        w = p * (y > z) + (1-p) * (y < z)
+    return z
 
 
 if len(sys.argv) < 3:
@@ -38,13 +38,18 @@ except OSError as e:
         raise
 
 
+# peak prominence threshold (for removing peaks and determining baseline)
+pk_prom_thres = 5.0
+# number of events to be plotted
+nev_plot = 10
+
 # sample data
 df = pd.read_csv(data_path, index_col=[0])
 # cuts on events
-dfc = df.loc[(df['Signal Sum'] < 1000) & (df['Fired Counts'] > 8)]
-# peak prominence threshold (for removing peaks and determining baseline)
-pk_prom_thres = 5.0
-for il in np.arange(3):
+# dfc = df.loc[(df['Signal Sum'] < 1000) & (df['Fired Counts'] > 8)]
+dfc = df.loc[df['LED'] > 0]
+
+for il in np.arange(nev_plot):
     if dfc.shape[0] <= il:
         continue
     data = dfc.iloc[il]
@@ -59,7 +64,9 @@ for il in np.arange(3):
         peaks, peak_props = signal.find_peaks(samples, prominence=pk_prom_thres)
         peak_intervals = pd.IntervalIndex.from_arrays(peak_props['left_bases'], peak_props['right_bases'])
         baseline = np.average(samples, weights=[~peak_intervals.contains(x).any() for x in np.arange(len(samples))])
+#        baseline = baseline_als(samples)
         ax.plot(np.arange(len(samples)), samples)
+#        ax.plot(np.arange(len(samples)), baseline, 'r--')
         ax.axhline(y=baseline, color='r')
 
     # fig.text(0.5, 0.04, 'Sample Number', ha='center', fontsize=24)

src/analyze.cpp

@@ -332,8 +332,8 @@ void process_data(const std::string &dpath, const std::string &opath, int nev, c
         }
 
         // further process event
-        // no reference time
-        if (event["Ref"].size() == 0)
+        // no reference time or triggered by LED
+        if ((event["Ref"].size() == 0) || (event["LED"].size() > 0))
             continue;
         auto ref_time = event["Ref"].front().time;
         // Calosum timing

src/extract_samples.cpp

@@ -192,7 +192,7 @@ void process_data(const std::string &dpath, const std::string &opath, int nev, c
     char csv_sep = ',';
 
     // csv file header
-    output << "Event Number" << csv_sep << "Fired Counts" << csv_sep <<"Signal Sum";
+    output << "Event Number" << csv_sep << "Fired Counts" << csv_sep <<"Signal Sum" << csv_sep << "LED";
     for (auto &sig : signals) {
         output << csv_sep << sig;
     }
@@ -261,8 +261,9 @@ void process_data(const std::string &dpath, const std::string &opath, int nev, c
             }
         }
 
-        if ((fire > 6) && has_sample(samples)) {
-            output << count << csv_sep << fire << csv_sep << sum;
+        bool led_trigger = (event["LED"].size() > 0);
+        if (led_trigger || ((fire > 6) && has_sample(samples))) {
+            output << count << csv_sep << fire << csv_sep << sum << csv_sep << led_trigger;
             for (auto &sig : signals) {
                 output << csv_sep << "\"[";
                 for (auto &spl : samples[sig]) {