|
10 | 10 | import update_path |
11 | 11 | from plt_gen import plt_gen |
12 | 12 | from get_params import get_params |
| 13 | +from astropy.stats import sigma_clipped_stats |
13 | 14 |
|
14 | 15 | import os |
15 | 16 |
|
16 | 17 | from pathlib import Path |
17 | 18 | import argparse |
18 | 19 |
|
19 | 20 | import matplotlib.pyplot as plt |
| 21 | +import matplotlib.ticker as ticker |
| 22 | + |
| 23 | +from scipy.optimize import curve_fit |
20 | 24 | import numpy as np |
| 25 | +import pandas as pd |
21 | 26 |
|
22 | 27 | from read import read_pflib_csv |
23 | 28 |
|
|
30 | 35 | parser.add_argument('-od', '--output_directory', type=Path, help='directory to which to print.') |
31 | 36 | parser.add_argument('-r', '--repository', type=bool, help='True if you would like to indicate that the input is a repository with csv files rather than a single csv file.') |
32 | 37 | plot_types = ['SCATTER', 'HEATMAP'] |
33 | | -plot_funcs = ['ADC-TIME', 'TOT-TIME', 'TOT', 'TOT-EFF', 'PARAMS', 'MULTI-CHANNEL'] |
| 38 | +plot_funcs = ['ADC-TIME', 'ADC-ALL-CHANNELS', 'TOT-TIME', 'TOT', 'TOT-EFF', 'PARAMS', 'MULTI-CHANNEL'] |
34 | 39 | parser.add_argument('-pt','--plot_type', choices=plot_types, default=plot_types[0], type=str, help=f'Plotting type. Options: {", ".join(plot_types)}') |
35 | 40 | parser.add_argument('-pf','--plot_function', choices=plot_funcs, default=plot_types[0], type=str, help=f'Plotting function. Options: {", ".join(plot_types)}') |
36 | 41 | parser.add_argument('-xl','--xlabel', default='time [ns]', type=str, help=f'What to label the x-axis with.') |
@@ -71,10 +76,197 @@ def set_xticks ( |
71 | 76 | xmin = 25*np.floor(xmin/25) |
72 | 77 | xmax = 25*np.ceil(xmax/25) |
73 | 78 | ax.set_xticks(ticks = np.arange(xmin, xmax+1, 25), minor=False) |
74 | | - ax.set_xticks(ticks = np.arange(xmin, xmax, 25/16), minor=True) |
| 79 | + ax.set_xticks(ticks = np.arange(xmin, xmax, 25/16), minor=True) |
| 80 | + |
| 81 | +#################### HELPER FUNCTIONS ######################## |
| 82 | + |
| 83 | +def linear(x, k, m): |
| 84 | + return k*x + m |
75 | 85 |
|
76 | 86 | #################### PLOTTING FUNCTIONS ######################## |
77 | 87 |
|
| 88 | +def adc_all_channels( |
| 89 | + samples, |
| 90 | + run_params, |
| 91 | + ax_holder, |
| 92 | + plot_params="all", # "all" or "one" |
| 93 | + param_index=6 # used only when plot_params="one" |
| 94 | +): |
| 95 | + """ |
| 96 | + Plot max ADC vs channels, compute RMS trends, and fit dependence on parameters. |
| 97 | + """ |
| 98 | + |
| 99 | + def compute_rms(ch_df, pedestal): |
| 100 | + """RMS for a single channel dataframe.""" |
| 101 | + return np.sqrt(((ch_df['adc'] - pedestal)**2).mean()) |
| 102 | + |
| 103 | + def compute_all_rms(df, pedestal): |
| 104 | + """Return {channel: rms}""" |
| 105 | + rms = df.groupby('channel').apply(lambda c: compute_rms(c, pedestal)) |
| 106 | + return rms.reset_index(name='rms') |
| 107 | + |
| 108 | + def conditional_legend(ax, max_labels=37, **kwargs): |
| 109 | + """Only draw legend if below threshold.""" |
| 110 | + handles, labels = ax.get_legend_handles_labels() |
| 111 | + if len(labels) < max_labels: |
| 112 | + ax.legend(**kwargs) |
| 113 | + |
| 114 | + def save_and_close(fig, fname): |
| 115 | + fig.savefig(fname, dpi=400) |
| 116 | + plt.close(fig) |
| 117 | + |
| 118 | + charges_group = samples.groupby("nr channels") |
| 119 | + fig_rms, ax_rms = plt.subplots(1, 1) |
| 120 | + ax_rms.set_xlabel('Channels activated') |
| 121 | + ax_rms.set_ylabel('Average RMS of non-activated channels') |
| 122 | + |
| 123 | + link = 0 |
| 124 | + central_val = 54 if link == 1 else 18 |
| 125 | + |
| 126 | + activated_channels_list = [] |
| 127 | + avg_rms_list = [] |
| 128 | + |
| 129 | + runs = len(charges_group) |
| 130 | + activated_channels_list = [[] for _ in range(runs - 1)] |
| 131 | + avg_rms_list = [[] for _ in range(runs - 1)] |
| 132 | + |
| 133 | + parameter_values = set() |
| 134 | + |
| 135 | + # Loop over charge groups |
| 136 | + for i, (charge_id, charge_df) in enumerate(charges_group): |
| 137 | + print(f'Running {i+1} out of {runs}') |
| 138 | + # Pedestal case |
| 139 | + if i == 0: |
| 140 | + fig, ax = plt.subplots() |
| 141 | + df = charge_df[charge_df['nr channels'] == 0] |
| 142 | + |
| 143 | + df = df[df['channel'] < 36] if link == 0 else df[df['channel'] >= 36] |
| 144 | + |
| 145 | + mean, med, std = sigma_clipped_stats(df['adc'], sigma=3) |
| 146 | + pedestal = mean |
| 147 | + |
| 148 | + ax.scatter(df['channel'], df['adc'], s=5, color='r', lw=1) |
| 149 | + ax.set_ylabel('ADC') |
| 150 | + ax.set_xlabel('Channel') |
| 151 | + |
| 152 | + save_and_close(fig, 'Pedestal.png') |
| 153 | + continue |
| 154 | + |
| 155 | + fig, ax1 = plt.subplots() |
| 156 | + fig_time, ax_time = plt.subplots() |
| 157 | + |
| 158 | + ax1.set_ylabel(r'$\Delta$ADC') |
| 159 | + ax1.set_xlabel('Channel') |
| 160 | + ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator()) |
| 161 | + |
| 162 | + ax_time.set_ylabel('ADC') |
| 163 | + ax_time.set_xlabel('Time') |
| 164 | + |
| 165 | + # Parameter grouping |
| 166 | + try: |
| 167 | + param_group, param_name = get_params(charge_df, 0) |
| 168 | + except Exception: |
| 169 | + param_group = [(None, charge_df)] |
| 170 | + |
| 171 | + cmap = plt.get_cmap('viridis') |
| 172 | + |
| 173 | + for j, (param_id, param_df) in enumerate(param_group): |
| 174 | + |
| 175 | + # --- Parameter selection --- |
| 176 | + if plot_params == "one": |
| 177 | + if j != param_index: |
| 178 | + continue |
| 179 | + |
| 180 | + try: |
| 181 | + key = param_name.split('.')[1] |
| 182 | + val = param_df[param_name].iloc[0] |
| 183 | + except Exception: |
| 184 | + key = "Channels flashed, voltage [mVpp]" |
| 185 | + val = (4, 2000) |
| 186 | + |
| 187 | + parameter_values.add(val) |
| 188 | + |
| 189 | + max_diff = ( |
| 190 | + param_df |
| 191 | + .groupby('channel')['adc'] |
| 192 | + .apply(lambda s: (s - pedestal).abs().max()) |
| 193 | + .reset_index(name='max_diff') |
| 194 | + ) |
| 195 | + |
| 196 | + # Get RMS |
| 197 | + rms_df = compute_all_rms(param_df, pedestal) |
| 198 | + rms_vals = rms_df['rms'].tolist() |
| 199 | + avg_rms = 0 |
| 200 | + count = 0 |
| 201 | + for k in range(central_val - 18, central_val + 18): |
| 202 | + if central_val - (i - 1) <= k <= central_val + (i - 1): |
| 203 | + continue |
| 204 | + avg_rms += rms_vals[k] |
| 205 | + count += 1 |
| 206 | + |
| 207 | + if count > 0: |
| 208 | + avg_rms /= count |
| 209 | + avg_rms_list[i - 1].append(avg_rms) |
| 210 | + activated_channels_list[i - 1].append(i) |
| 211 | + |
| 212 | + merged = max_diff.merge(rms_df, on='channel') |
| 213 | + ax1.set_ylim(-10, merged['max_diff'].max() + 10) |
| 214 | + ax1.scatter( |
| 215 | + merged['channel'], merged['max_diff'], |
| 216 | + label=f'{key} = {val}', |
| 217 | + s=5, color=cmap(j / len(param_group)), lw=1 |
| 218 | + ) |
| 219 | + if link == 0: |
| 220 | + link_df = param_df[param_df['channel'] < 36] |
| 221 | + elif link == 1: |
| 222 | + link_df = param_df[param_df['channel'] >= 36] |
| 223 | + for k, (ch_id, ch_df) in enumerate(link_df.groupby('channel')): |
| 224 | + ax_time.scatter( |
| 225 | + ch_df['time'], ch_df['adc'], |
| 226 | + label=f'ch{ch_id}', |
| 227 | + s=5, color=plt.get_cmap('tab20')(k / 20) |
| 228 | + ) |
| 229 | + |
| 230 | + # Pedestal and link lines |
| 231 | + ax1.axhline(y=0, color='k', linestyle='--', linewidth=0.8) |
| 232 | + ax1.axvline(x=35.5, color='r', linestyle='--', alpha=0.5) |
| 233 | + |
| 234 | + conditional_legend(ax1, fontsize=8) |
| 235 | + conditional_legend(ax_time, fontsize=4, ncols=3) |
| 236 | + |
| 237 | + save_and_close(fig, f'adc_channels_{i}.png') |
| 238 | + save_and_close(fig_time, f'adc_time_{i}.png') |
| 239 | + |
| 240 | + # slope/intercept plots |
| 241 | + activated_T = list(map(list, zip(*activated_channels_list))) |
| 242 | + rms_T = list(map(list, zip(*avg_rms_list))) |
| 243 | + parameter_values = sorted(parameter_values) |
| 244 | + |
| 245 | + if len(activated_T) > 1: |
| 246 | + slopes = [] |
| 247 | + intercepts = [] |
| 248 | + |
| 249 | + for xvals, yvals in zip(activated_T, rms_T): |
| 250 | + popt, _ = curve_fit(linear, xvals, yvals) |
| 251 | + slopes.append(popt[0]) |
| 252 | + intercepts.append(popt[1]) |
| 253 | + |
| 254 | + ax_rms.scatter(xvals, yvals, s=8) |
| 255 | + xfit = np.linspace(min(xvals), max(xvals), 100) |
| 256 | + ax_rms.plot(xfit, linear(xfit, *popt), |
| 257 | + label=f'Fit: y={popt[0]:.3f}x + {popt[1]:.3f}, CALIB={parameter_values[len(slopes)-1]}') |
| 258 | + |
| 259 | + conditional_legend(ax_rms, fontsize=8) |
| 260 | + save_and_close(fig_rms, 'avg_rms.png') |
| 261 | + |
| 262 | + fig_par, ax_par = plt.subplots() |
| 263 | + ax_par.scatter(parameter_values, slopes, label='slopes') |
| 264 | + ax_par.scatter(parameter_values, intercepts, label='intercepts') |
| 265 | + ax_par.legend() |
| 266 | + ax_par.set_xlabel('CALIB') |
| 267 | + ax_par.set_title('Slope and intercepts from fits') |
| 268 | + save_and_close(fig_par, 'slope_intercept.png') |
| 269 | + |
78 | 270 | def time( |
79 | 271 | samples, |
80 | 272 | run_params, |
@@ -374,6 +566,10 @@ def multiparams( |
374 | 566 |
|
375 | 567 | ############################## MAIN ################################### |
376 | 568 |
|
| 569 | +if args.plot_function == 'ADC-ALL-CHANNELS': |
| 570 | + plt_gen(adc_all_channels, samples, run_params, args.output, |
| 571 | + xlabel = args.xlabel, ylabel = args.ylabel) |
| 572 | + |
377 | 573 | if args.plot_function == 'ADC-TIME' and args.plot_type == 'SCATTER': |
378 | 574 | plt_gen(partial(time, yval = 'adc'), samples, run_params, args.output, |
379 | 575 | xlabel = args.xlabel, ylabel = args.ylabel) |
|
0 commit comments