From f5f8a6ab1c047672422e5fdeb6012f5f020f92e5 Mon Sep 17 00:00:00 2001 From: YannickDieter Date: Fri, 17 Nov 2017 16:11:28 +0100 Subject: [PATCH 1/4] ENH: add in-pixel-hit-distribution for CS and eff. pitch calculation --- testbeam_analysis/result_analysis.py | 270 ++++++++++++++++++++++ testbeam_analysis/tools/analysis_utils.py | 39 ++++ testbeam_analysis/tools/plot_utils.py | 204 ++++++++++++++++ 3 files changed, 513 insertions(+) diff --git a/testbeam_analysis/result_analysis.py b/testbeam_analysis/result_analysis.py index 90e1ae8a..ec088d64 100644 --- a/testbeam_analysis/result_analysis.py +++ b/testbeam_analysis/result_analysis.py @@ -1299,3 +1299,273 @@ def histogram_track_angle(input_tracks_file, input_alignment_file=None, output_t if plot: plot_utils.plot_track_angle(input_track_angle_file=output_track_angle_file, output_pdf_file=None, dut_names=dut_names) + + +def calculate_in_pixel_hit_distribution(input_tracks_file, input_alignment_file, use_duts, pixel_size, n_bins, output_pdf_file=None, plot=True, force_prealignment=False, output_in_pixel_dist_file=None, chunk_size=1000000): + '''Takes the tracks and calculates in-pixel-hit-distributions for cluster sizes between 1 and 4, additionally the effective CS-1-Pitch is calculated. + The effective CS-1-pitch 'p_eff' gives the area of one pixel in which CS1 is formed depending on the sensor threshold. It can be used to estimate the + intrinsic resolution: sigma_int = p_eff / sqrt(12) + + Parameters + ---------- + input_tracks_file : string + Filename of the input tracks file. + input_alignment_file : string + Filename of the input alignment file. + use_duts : iterable + Calculate in-pixel-hit-distribution for selected DUTs. If None, all duts are selected. + pixel_size : iterable of tuples + One tuple per DUT describing the pixel dimension in um in column, row direction + e.g. for 2 DUTs: pixel_size = [(250, 50), (250, 50)] + n_bins : iterable + Number of bins for the histogram in x and y direction. + output_in_pixel_dist_file : string + Filename of the output in-pixel-hit-distribution file. + plot : bool + If True, create output plots. If false, only effective CS-1-Pitch will be calculated + force_prealignment : boolean + Take the prealignment, although if a coarse alignment is availale. + chunk_size : uint + Chunk size of the data when reading from file. + ''' + logging.info('=== Calculating In-Pixel-Hit-Distribution for CS 1 - 4 ===') + + use_prealignment = True if force_prealignment else False + + def calculate_effective_cs1_pitch(cs_hist, pitch): + '''Calculates the effectice CS-1 pitch p_eff using the CS distribution: p_eff = sqrt(pitch^2 * CS1_percentage) + ''' + perc_cluster_size_1 = 100. * ((cs_hist[0].astype(np.float) / np.sum(cs_hist).astype(np.float))) + return np.sqrt(pitch * pitch * perc_cluster_size_1 / 100.) + + def normalize_distributions(hists, fit_results): + '''Normalize CS distributions in order to compare them. + ''' + hist_normed = [hists[0] - fit_results[0][3] + fit_results[2][0] + fit_results[2][3], # shift CS 1 distribution (x-axis projection) such that backgrounds of CS 2 and CS 1 is the same + hists[1], # CS 2 distribution (x-axis projection) stays the same + hists[2] - fit_results[1][3] + fit_results[3][0] + fit_results[3][3], # shift CS 1 distribution (y-axis projection) such that backgrounds of CS 2 and CS 1 is the same + hists[3]] # CS 2 distribution (y-axis projection) stays the same + return hist_normed + + # read alignment for converting the track intersections with the respective dut (global reference system) to local track intersections + with tb.open_file(input_alignment_file, mode="r") as in_file_h5: # Open file with alignment data + if use_prealignment: + logging.info('Use pre-alignment data') + prealignment = in_file_h5.root.PreAlignment[:] + n_duts = prealignment.shape[0] + else: + logging.info('Use alignment data') + alignment = in_file_h5.root.Alignment[:] + n_duts = alignment.shape[0] + + use_duts = use_duts if use_duts is not None else range(n_duts) # standard setting: fit tracks for all DUTs + + if plot: + output_pdf = PdfPages(output_pdf_file + '.pdf') + else: + output_pdf = None + + n_bins = [n_bins, n_bins] if not isinstance(n_bins, Iterable) else n_bins + + effective_pitches = [] + + with tb.open_file(input_tracks_file, 'r') as in_file_h5: + for actual_dut in use_duts: + node = in_file_h5.get_node(in_file_h5.root, 'Tracks_DUT_%d' % actual_dut) + + actual_pixel_size = pixel_size[actual_dut] + + initialize = True + for tracks_chunk, _ in analysis_utils.data_aligned_at_events(node, chunk_size=chunk_size): # read track file in chunks + print tracks_chunk.dtype + # select only valid track intersection + tracks_chunk = tracks_chunk[~np.isnan(tracks_chunk[:]['offset_0'])] + + # select different CS in order to histogram them and to deduce from the CS ratio the effective CS-1-pitch + n_hits = tracks_chunk[~np.isnan(tracks_chunk['n_hits_dut_%d' % actual_dut])]['n_hits_dut_%d' % actual_dut] + + # read track intersections wiht actual dut for cluster sizes betwenn 1 and 4 + intersection_x, intersection_y, intersection_z = tracks_chunk[:]['offset_0'], tracks_chunk[:]['offset_1'], tracks_chunk[:]['offset_2'] + intersection_x_cs_1 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] + intersection_y_cs_1 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] + intersection_z_cs_1 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] + intersection_x_cs_2 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] + intersection_y_cs_2 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] + intersection_z_cs_2 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] + intersection_x_cs_3 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] + intersection_y_cs_3 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] + intersection_z_cs_3 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] + intersection_x_cs_4 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] + intersection_y_cs_4 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] + intersection_z_cs_4 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] + + # transoform to local coordinate system + if use_prealignment: + intersection_x_cs_1_local, intersection_y_cs_1_local, intersection_z_cs_1_local = geometry_utils.apply_alignment(intersection_x_cs_1, intersection_y_cs_1, intersection_z_cs_1, + dut_index=actual_dut, + prealignment=prealignment, + inverse=True) + + intersection_x_cs_2_local, intersection_y_cs_2_local, intersection_z_cs_2_local = geometry_utils.apply_alignment(intersection_x_cs_2, intersection_y_cs_2, intersection_z_cs_2, + dut_index=actual_dut, + prealignment=prealignment, + inverse=True) + + intersection_x_cs_3_local, intersection_y_cs_3_local, intersection_z_cs_3_local = geometry_utils.apply_alignment(intersection_x_cs_3, intersection_y_cs_3, intersection_z_cs_3, + dut_index=actual_dut, + prealignment=prealignment, + inverse=True) + + intersection_x_cs_4_local, intersection_y_cs_4_local, intersection_z_cs_4_local = geometry_utils.apply_alignment(intersection_x_cs_4, intersection_y_cs_4, intersection_z_cs_4, + dut_index=actual_dut, + prealignment=prealignment, + inverse=True) + + else: + intersection_x_cs_1_local, intersection_y_cs_1_local, intersection_z_cs_1_local = geometry_utils.apply_alignment(intersection_x_cs_1, intersection_y_cs_1, intersection_z_cs_1, + dut_index=actual_dut, + alignment=alignment, + inverse=True) + + intersection_x_cs_2_local, intersection_y_cs_2_local, intersection_z_cs_2_local = geometry_utils.apply_alignment(intersection_x_cs_2, intersection_y_cs_2, intersection_z_cs_2, + dut_index=actual_dut, + alignment=alignment, + inverse=True) + + intersection_x_cs_3_local, intersection_y_cs_3_local, intersection_z_cs_3_local = geometry_utils.apply_alignment(intersection_x_cs_3, intersection_y_cs_3, intersection_z_cs_3, + dut_index=actual_dut, + alignment=alignment, + inverse=True) + + intersection_x_cs_4_local, intersection_y_cs_4_local, intersection_z_cs_4_local = geometry_utils.apply_alignment(intersection_x_cs_4, intersection_y_cs_4, intersection_z_cs_4, + dut_index=actual_dut, + alignment=alignment, + inverse=True) + + if not np.allclose(intersection_z_cs_4_local[np.isfinite(intersection_z_cs_4_local)], 0.0): + print intersection_z_cs_4_local + raise RuntimeError("Transformation into local coordinate system gives z != 0") + + # project track intersections onto one pixel in order to increase statistics + projection_x_1 = np.mod(intersection_x_cs_1_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_1_local))) + projection_y_1 = np.mod(intersection_y_cs_1_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_1_local))) + projection_x_2 = np.mod(intersection_x_cs_2_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_2_local))) + projection_y_2 = np.mod(intersection_y_cs_2_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_2_local))) + projection_x_3 = np.mod(intersection_x_cs_3_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_3_local))) + projection_y_3 = np.mod(intersection_y_cs_3_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_3_local))) + projection_x_4 = np.mod(intersection_x_cs_4_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_4_local))) + projection_y_4 = np.mod(intersection_y_cs_4_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_4_local))) + + # histogram intersections and CS (for calculation of effective pitch) + if initialize: + projection_x_1_hist, edges_x = np.histogram(projection_x_1, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projection_y_1_hist, edges_y = np.histogram(projection_y_1, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + projection_x_2_hist, _ = np.histogram(projection_x_2, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projection_y_2_hist, _ = np.histogram(projection_y_2, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + projection_x_3_hist, _ = np.histogram(projection_x_3, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projection_y_3_hist, _ = np.histogram(projection_y_3, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + projection_x_4_hist, _ = np.histogram(projection_x_4, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projection_y_4_hist, _ = np.histogram(projection_y_4, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + + projection_1_2d_hist, _, _ = np.histogram2d(x=projection_x_1, y=projection_y_1, bins=[edges_x, edges_y]) + projection_2_2d_hist, _, _ = np.histogram2d(x=projection_x_2, y=projection_y_2, bins=[edges_x, edges_y]) + projection_3_2d_hist, _, _ = np.histogram2d(x=projection_x_3, y=projection_y_3, bins=[edges_x, edges_y]) + projection_4_2d_hist, _, _ = np.histogram2d(x=projection_x_4, y=projection_y_4, bins=[edges_x, edges_y]) + + cs_edges = np.arange(1.0 - 0.5, 1000 + 0.5, 1) + cs_hist, cs_edges = np.histogram(n_hits, bins=cs_edges, density=False) + initialize = False + else: # if already read first chunk, add histograms up + projection_x_1_hist += np.histogram(projection_x_1, bins=edges_x)[0] + projection_y_1_hist += np.histogram(projection_y_1, bins=edges_y)[0] + projection_x_2_hist += np.histogram(projection_x_2, bins=edges_x)[0] + projection_y_2_hist += np.histogram(projection_y_2, bins=edges_y)[0] + projection_x_3_hist += np.histogram(projection_x_3, bins=edges_x)[0] + projection_y_3_hist += np.histogram(projection_y_3, bins=edges_y)[0] + projection_x_4_hist += np.histogram(projection_x_4, bins=edges_x)[0] + projection_y_4_hist += np.histogram(projection_y_4, bins=edges_y)[0] + + projection_1_2d_hist += np.histogram2d(x=projection_x_1, y=projection_y_1, bins=[edges_x, edges_y])[0] + projection_2_2d_hist += np.histogram2d(x=projection_x_2, y=projection_y_2, bins=[edges_x, edges_y])[0] + projection_3_2d_hist += np.histogram2d(x=projection_x_3, y=projection_y_3, bins=[edges_x, edges_y])[0] + projection_4_2d_hist += np.histogram2d(x=projection_x_4, y=projection_y_4, bins=[edges_x, edges_y])[0] + + cs_hist += np.histogram(n_hits, bins=cs_edges, density=False)[0] + + # calculate effective pitch from cluster size ratio + effective_pitch = calculate_effective_cs1_pitch(cs_hist=cs_hist, pitch=actual_pixel_size[0]) + effective_pitches.append(effective_pitch) + + if plot: + # plot histograms for cluster sizes between 1 and 4 + plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_1_hist, + y_intersections_hist=projection_y_1_hist, + intersections_2d_hist=projection_1_2d_hist, + pixel_pitch=actual_pixel_size, + output_pdf=output_pdf, + bins=[edges_x, edges_y], + plot_title=('In-pixel Cluster Size 1 Hit Distribution for DUT%d' % actual_dut)) + plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_2_hist, + y_intersections_hist=projection_y_2_hist, + intersections_2d_hist=projection_2_2d_hist, + pixel_pitch=actual_pixel_size, + output_pdf=output_pdf, + bins=[edges_x, edges_y], + plot_title=('In-pixel Cluster Size 2 Hit Distribution for DUT%d' % actual_dut)) + plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_3_hist, + y_intersections_hist=projection_y_3_hist, + intersections_2d_hist=projection_3_2d_hist, + pixel_pitch=actual_pixel_size, + output_pdf=output_pdf, + bins=[edges_x, edges_y], + plot_title=('In-pixel Cluster Size 3 Hit Distribution for DUT%d' % actual_dut)) + plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_4_hist, + y_intersections_hist=projection_y_4_hist, + intersections_2d_hist=projection_4_2d_hist, + pixel_pitch=actual_pixel_size, + output_pdf=output_pdf, + bins=[edges_x, edges_y], + plot_title=('In-pixel Cluster Size 4 Hit Distribution for DUT%d' % actual_dut)) + + # take slice out of middle + slice_width = 5 # set withd of slice to 5 um + slice_start_x = (actual_pixel_size[0] - slice_width) / 2.0 + slice_stop_x = (actual_pixel_size[0] + slice_width) / 2.0 + slice_start_y = (actual_pixel_size[1] - slice_width) / 2.0 + slice_stop_y = (actual_pixel_size[1] + slice_width) / 2.0 + start_x = int(slice_start_x / (actual_pixel_size[0] / n_bins[0])) + stop_x = int(slice_stop_x / (actual_pixel_size[0] / n_bins[0])) + start_y = int(slice_start_y / (actual_pixel_size[0] / n_bins[1])) + stop_y = int(slice_stop_y / (actual_pixel_size[0] / n_bins[1])) + sliced_hist_y_1 = np.sum(projection_1_2d_hist[start_x:stop_x, ], axis=0) + sliced_hist_x_1 = np.sum(projection_1_2d_hist[:, start_y:stop_y], axis=1) + sliced_hist_y_2 = np.sum(projection_2_2d_hist[start_x:stop_x, ], axis=0) + sliced_hist_x_2 = np.sum(projection_2_2d_hist[:, start_y:stop_y], axis=1) + + # fit projections onto x and y-axis of CS 1 and CS 2 distributions + fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced = analysis_utils.fit_in_pixel_hit_hist(cs1_hist_projected=[sliced_hist_x_1, sliced_hist_y_1], + cs2_hist_projected=[sliced_hist_x_2, sliced_hist_y_2], + bin_center_x=(edges_x[1:] + edges_x[:-1]) / 2.0, + bin_center_y=(edges_y[1:] + edges_y[:-1]) / 2.0) + + # deduce normalization (shift) from fit, norm such that background of CS 1 distribution is the same as background of CS 2 distribution + x_intersections_hist_cs_1_normed, x_intersections_hist_cs_2_normed, y_intersections_hist_cs_1_normed, y_intersections_hist_cs_2_normed = normalize_distributions(hists=[sliced_hist_x_1, sliced_hist_x_2, sliced_hist_y_1, sliced_hist_y_2], + fit_results=[fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced]) + if plot: + # plot effective CS-1 pitch from ratio of CS distribution + plot_utils.plot_in_pixel_hit_hist_with_eff_pitch(x_intersections_hist_cs_1=x_intersections_hist_cs_1_normed, + x_intersections_hist_cs_2=x_intersections_hist_cs_2_normed, + y_intersections_hist_cs_1=y_intersections_hist_cs_1_normed, + y_intersections_hist_cs_2=y_intersections_hist_cs_2_normed, + intersections_2d_hist=projection_1_2d_hist - projection_2_2d_hist, + fit_results=[fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced], + pixel_pitch=actual_pixel_size, + bins=[edges_x, edges_y], + effective_pitch=effective_pitch, + output_pdf=output_pdf, + plot_title=('Effective Cluster Size 1 Pitch for DUT%d' % actual_dut)) + + if output_pdf is not None: + output_pdf.close() + + return effective_pitches diff --git a/testbeam_analysis/tools/analysis_utils.py b/testbeam_analysis/tools/analysis_utils.py index 539705ee..ce67bf27 100644 --- a/testbeam_analysis/tools/analysis_utils.py +++ b/testbeam_analysis/tools/analysis_utils.py @@ -924,6 +924,45 @@ def fit_residuals_vs_position(hist, xedges, yedges, xlabel="", ylabel="", title= return fit, cov +def fit_in_pixel_hit_hist(cs1_hist_projected, cs2_hist_projected, bin_center_x, bin_center_y): + # fit distributions with gauss curve + fit_cs_1_x_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, + bin_center_x, + cs1_hist_projected[0], + p0=(np.max(cs1_hist_projected[0]) - np.min(cs1_hist_projected[0]), + np.median(bin_center_x), + np.std(np.repeat((bin_center_x).astype(np.int64), (cs1_hist_projected[0]).astype(np.int64))), + np.min(cs1_hist_projected[0])), + sigma=np.sqrt(cs1_hist_projected[0])) + fit_cs_1_y_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, + bin_center_y, + cs1_hist_projected[1], + p0=(np.max(cs1_hist_projected[1]) - np.min(cs1_hist_projected[1]), + np.median(bin_center_y), + np.std(np.repeat((bin_center_y).astype(np.int64), (cs1_hist_projected[1]).astype(np.int64))), + np.min(cs1_hist_projected[1])), + sigma=np.sqrt(cs1_hist_projected[1])) + + fit_cs_2_x_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, + bin_center_x, + cs2_hist_projected[0], + p0=(np.max(cs2_hist_projected[0]) - np.min(cs2_hist_projected[0]), + np.median(bin_center_x), + np.std(np.repeat((bin_center_x).astype(np.int64), (cs2_hist_projected[0]).astype(np.int64))), + np.min(cs2_hist_projected[0])), + sigma=np.sqrt(cs2_hist_projected[0])) + fit_cs_2_y_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, + bin_center_y, + cs2_hist_projected[1], + p0=(np.max(cs2_hist_projected[1]) - np.min(cs2_hist_projected[1]), + np.median(bin_center_y), + np.std(np.repeat((bin_center_y).astype(np.int64), (cs2_hist_projected[1]).astype(np.int64))), + np.min(cs2_hist_projected[1])), + sigma=np.sqrt(cs2_hist_projected[1])) + + return fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced + + def hough_transform(img, theta_res=1.0, rho_res=1.0, return_edges=False): thetas = np.linspace(-90.0, 0.0, np.ceil(90.0/theta_res) + 1) thetas = np.concatenate((thetas, -thetas[len(thetas)-2::-1])) diff --git a/testbeam_analysis/tools/plot_utils.py b/testbeam_analysis/tools/plot_utils.py index 0d721b0e..84d3da1c 100644 --- a/testbeam_analysis/tools/plot_utils.py +++ b/testbeam_analysis/tools/plot_utils.py @@ -9,6 +9,8 @@ import numpy as np import tables as tb import matplotlib as mpl +import matplotlib.gridspec as gridspec +from matplotlib import ticker from matplotlib.backends.backend_pdf import PdfPages from matplotlib.backends.backend_agg import FigureCanvas import matplotlib.pyplot as plt @@ -18,12 +20,24 @@ from mpl_toolkits.mplot3d import Axes3D # needed for 3d plotting although it is shown as not used from matplotlib.widgets import Slider, Button from scipy.optimize import curve_fit +from scipy.optimize import fsolve import testbeam_analysis.tools.analysis_utils warnings.filterwarnings("ignore", category=UserWarning, module="matplotlib") # Plot backend error not important +# class to define the middle of a colorbar +class MidpointNormalize(colors.Normalize): + def __init__(self, vmin=None, vmax=None, midpoint=None, clip=False): + self.midpoint = midpoint + colors.Normalize.__init__(self, vmin, vmax, clip) + + def __call__(self, value, clip=None): + x, y = [self.vmin, self.midpoint, self.vmax], [0, 0.5, 1] + return np.ma.masked_array(np.interp(value, x, y)) + + def plot_2d_pixel_hist(fig, ax, hist2d, plot_range, title=None, x_axis_title=None, y_axis_title=None, z_min=0, z_max=None): extent = [0.5, plot_range[0] + .5, plot_range[1] + .5, 0.5] if z_max is None: @@ -1333,3 +1347,193 @@ def plot_track_angle(input_track_angle_file, output_pdf_file=None, dut_names=Non ax.grid() ax.set_xlim(min(edges), max(edges)) output_pdf.savefig(fig) + + +def plot_in_pixel_hit_hist(x_intersections_hist, y_intersections_hist, intersections_2d_hist, pixel_pitch, bins, output_pdf=None, plot_title=None, z_max=None): + '''Plotting in-pixel-hit-distributions for cluster sizes between 1 and 4. + ''' + logging.info('Plotting In-Pixel-Hit-Distribution for CS 1 - 4') + if not output_pdf: + return + if z_max is None: + if intersections_2d_hist.all() is np.ma.masked: # check if masked array is fully masked + z_max = 1 + else: + z_max = ceil(intersections_2d_hist.max()) + + bin_center_x = (bins[0][1:] + bins[0][:-1]) / 2.0 # center of bin + bin_center_y = (bins[1][1:] + bins[1][:-1]) / 2.0 # center of bin + + fig = Figure(figsize=(12., 12.)) + _ = FigureCanvas(fig) + + # colors + cmap = plt.cm.get_cmap('viridis') + plot_color = mpl.colors.to_hex(cmap(.2), keep_alpha=True) + + # setup plot + gs1 = gridspec.GridSpec(2, 2, width_ratios=[3, 1], height_ratios=[1, 3]) + gs1.update(wspace=0.05, hspace=0.05) + axHistx = fig.add_subplot(gs1[:-1, :-1]) + axHisty = fig.add_subplot(gs1[-1, -1]) + axHeatmap = fig.add_subplot(gs1[-1, :-1]) + + axHisty.get_yaxis().set_ticklabels([]) + axHistx.get_xaxis().set_ticklabels([]) + axHistx.set_ylabel('#', fontsize=14) + axHisty.set_xlabel('#', fontsize=14) + axHeatmap.set_xlabel('$x_\mathrm{track}$ / $\mathrm{\mu}$m', fontsize=14) + axHeatmap.set_ylabel('$y_\mathrm{track}$ / $\mathrm{\mu}$m', fontsize=14) + axHeatmap.set_aspect('equal') + _ = [tick.label.set_fontsize(14) for tick in axHistx.yaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHisty.xaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHeatmap.yaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHeatmap.xaxis.get_major_ticks()] + nticks = plt.MaxNLocator(4) + formatter = ticker.ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((-1, 1)) + axHistx.yaxis.set_major_formatter(formatter) + axHistx.yaxis.set_major_locator(nticks) + axHisty.xaxis.set_major_formatter(formatter) + axHisty.xaxis.set_major_locator(nticks) + axHeatmap.legend(loc='upper right') + if plot_title is not None: + axHistx.set_title(plot_title) + + # plot in pixel hit distribution and projections on top of axes + heatmap_cs = axHeatmap.imshow(intersections_2d_hist, origin='lower', interpolation='none', aspect="auto", extent=[0.0, bins[0][-1], 0.0, bins[1][-1]], cmap=cmap, clim=(0, z_max)) + axHistx.bar(bin_center_x, x_intersections_hist, color=plot_color, width=np.mean(np.diff(bins[0]))) + axHisty.barh(bin_center_y, y_intersections_hist, color=plot_color, height=np.mean(np.diff(bins[1]))) + + # colorbar + cbar = fig.colorbar(heatmap_cs, ax=axHisty) + cbar.ax.tick_params(labelsize=14) + + # adjust plot limits + axHistx.set_xlim(0., pixel_pitch[0]) + axHisty.set_ylim(0., pixel_pitch[1]) + axHistx.set_ylim(bottom=0.8 * x_intersections_hist.min()) + axHisty.set_xlim(left=0.8 * y_intersections_hist.min()) + + output_pdf.savefig(fig) + + +def plot_in_pixel_hit_hist_with_eff_pitch(x_intersections_hist_cs_1, x_intersections_hist_cs_2, y_intersections_hist_cs_1, y_intersections_hist_cs_2, intersections_2d_hist, fit_results, pixel_pitch, bins, effective_pitch, output_pdf=None, plot_title=None): + '''Plotting in-pixel-hit-distributions for cluster sizes between 1 and 4. + ''' + logging.info('Plotting Effective CS-1-Pitch') + if not output_pdf: + return + + bin_center_x = (bins[0][1:] + bins[0][:-1]) / 2.0 # center of bin + bin_center_y = (bins[1][1:] + bins[1][:-1]) / 2.0 # center of bin + + fig = Figure(figsize=(12., 12.)) + _ = FigureCanvas(fig) + + # colors + cmap = plt.cm.get_cmap('RdBu') + plot_color_1 = mpl.colors.to_hex(cmap(.9), keep_alpha=True) + plot_color_2 = mpl.colors.to_hex(cmap(.1), keep_alpha=True) + + # setup plot + gs1 = gridspec.GridSpec(2, 2, width_ratios=[3, 1], height_ratios=[1, 3]) + gs1.update(wspace=0.05, hspace=0.05) + axHistx = fig.add_subplot(gs1[:-1, :-1]) + axHisty = fig.add_subplot(gs1[-1, -1]) + axHeatmap = fig.add_subplot(gs1[-1, :-1]) + + axHisty.get_yaxis().set_ticklabels([]) + axHistx.get_xaxis().set_ticklabels([]) + axHistx.set_ylabel('#', fontsize=14) + axHisty.set_xlabel('#', fontsize=14) + axHeatmap.set_xlabel('$x_\mathrm{track}$ / $\mathrm{\mu}$m', fontsize=14) + axHeatmap.set_ylabel('$y_\mathrm{track}$ / $\mathrm{\mu}$m', fontsize=14) + axHeatmap.set_aspect('equal') + _ = [tick.label.set_fontsize(14) for tick in axHistx.yaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHisty.xaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHeatmap.yaxis.get_major_ticks()] + _ = [tick.label.set_fontsize(14) for tick in axHeatmap.xaxis.get_major_ticks()] + nticks = plt.MaxNLocator(4) + formatter = ticker.ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((-1, 1)) + axHistx.yaxis.set_major_formatter(formatter) + axHistx.yaxis.set_major_locator(nticks) + axHisty.xaxis.set_major_formatter(formatter) + axHisty.xaxis.set_major_locator(nticks) + + if plot_title is not None: + axHistx.set_title(plot_title) + + x_gauss_x = np.arange(0.0, pixel_pitch[0], step=0.001) + x_gauss_y = np.arange(0.0, pixel_pitch[1], step=0.001) + + # unpack fit results + fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced = fit_results + + plot_fit_x_1_normed = testbeam_analysis.tools.analysis_utils.gauss_offset(x_gauss_x, *fit_cs_1_x_sliced) - fit_cs_1_x_sliced[3] + fit_cs_2_x_sliced[0] + fit_cs_2_x_sliced[3] + plot_fit_x_2_normed = testbeam_analysis.tools.analysis_utils.gauss_offset(x_gauss_x, *fit_cs_2_x_sliced) + plot_fit_y_1_normed = testbeam_analysis.tools.analysis_utils.gauss_offset(x_gauss_y, *fit_cs_1_y_sliced) - fit_cs_1_y_sliced[3] + fit_cs_2_y_sliced[0] + fit_cs_2_y_sliced[3] + plot_fit_y_2_normed = testbeam_analysis.tools.analysis_utils.gauss_offset(x_gauss_y, *fit_cs_2_y_sliced) + + # plot in pixel hit distribution and projections on top of axes and the respective fits + axHistx.bar(bin_center_x, x_intersections_hist_cs_1, width=np.mean(np.diff(bins[0])), color=plot_color_1) + axHistx.bar(bin_center_x, x_intersections_hist_cs_2, width=np.mean(np.diff(bins[0])), color=plot_color_2) + axHisty.barh(bin_center_y, y_intersections_hist_cs_1, height=np.mean(np.diff(bins[1])), color=plot_color_1) + axHisty.barh(bin_center_y, y_intersections_hist_cs_2, height=np.mean(np.diff(bins[1])), color=plot_color_2) + axHistx.plot(x_gauss_x, + plot_fit_x_1_normed, + color=mpl.colors.to_hex(cmap(1.), keep_alpha=True), + label='Cluster Size 1 Hit Distribution') + axHistx.plot(x_gauss_x, + plot_fit_x_2_normed, + color=mpl.colors.to_hex(cmap(0.), keep_alpha=True), + label='Cluster Size 2 Hit Distribution') + axHisty.plot(plot_fit_y_1_normed, + x_gauss_y, + color=mpl.colors.to_hex(cmap(1.), keep_alpha=True), + label='Cluster Size 1 Hit Distribution') + axHisty.plot(plot_fit_y_2_normed, + x_gauss_y, + color=mpl.colors.to_hex(cmap(0.), keep_alpha=True), + label='Cluster Size 2 Hit Distribution') + heatmap_cs = axHeatmap.imshow(intersections_2d_hist.T, origin='lower', interpolation='none', aspect="auto", extent=[0.0, bins[0][-1], 0.0, bins[1][-1]], cmap=cmap, norm=MidpointNormalize(midpoint=0., vmin=np.min(intersections_2d_hist), vmax=np.max(intersections_2d_hist))) + + # calculate intersections + def difference_x_sliced(x): + return testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_1_x_sliced) - fit_cs_1_x_sliced[3] + fit_cs_2_x_sliced[0] + fit_cs_2_x_sliced[3] - testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_2_x_sliced) + + def difference_y_sliced(x): + return testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_1_y_sliced) - fit_cs_1_y_sliced[3] + fit_cs_2_y_sliced[0] + fit_cs_2_y_sliced[3] - testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_2_y_sliced) + + roots_x_sliced = fsolve(difference_x_sliced, [(pixel_pitch[0] - effective_pitch) / 2.0, pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.0)]) + roots_y_sliced = fsolve(difference_y_sliced, [(pixel_pitch[1] - effective_pitch) / 2.0, pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.0)]) + + # draw intersections of CS 1 and CS 2 distribution + for i in range(len(roots_x_sliced)): + axHistx.axvline(roots_x_sliced[i], linestyle='--', color='dimgrey', label='Intersections' if i == 0 else None) + axHisty.axhline(roots_y_sliced[i], linestyle='--', color='dimgrey', label='Intersections' if i == 0 else None) + # draw effective cluster size area from overlap of CS1 and CS2 + axHeatmap.plot([roots_x_sliced[0], roots_x_sliced[0], roots_x_sliced[1], roots_x_sliced[1], roots_x_sliced[0]], + [roots_y_sliced[0], roots_y_sliced[1], roots_y_sliced[1], roots_y_sliced[0], roots_y_sliced[0]], + '-', color='black', label='Effective CS-1-Pitch from Transition of Cluster Size Hit Distributions') + # draw effective cluster size area from cluster size ratio + axHeatmap.plot([(pixel_pitch[0] - effective_pitch) / 2., (pixel_pitch[0] - effective_pitch) / 2., pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.), pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.), (pixel_pitch[0] - effective_pitch) / 2.], + [(pixel_pitch[1] - effective_pitch) / 2., pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.), pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.), (pixel_pitch[1] - effective_pitch) / 2., (pixel_pitch[1] - effective_pitch) / 2], + '--', color='black', label='Effective CS-1-Pitch from Cluster Size Ratio of Cluster Size Distribution') + + # colorbar + cbar = fig.colorbar(heatmap_cs, ax=axHisty) + cbar.ax.tick_params(labelsize=14) + + # adjust plot limits + axHistx.set_xlim(0., pixel_pitch[0]) + axHisty.set_ylim(0., pixel_pitch[1]) + axHistx.set_ylim(bottom=0.8 * x_intersections_hist_cs_2.min()) + axHisty.set_xlim(left=0.8 * y_intersections_hist_cs_2.min()) + # set legend + axHeatmap.legend(loc='upper right') + + output_pdf.savefig(fig) From 884ed022232d0892e017181642e864c8705ca975 Mon Sep 17 00:00:00 2001 From: YannickDieter Date: Thu, 7 Dec 2017 12:24:42 +0100 Subject: [PATCH 2/4] MAINT: more compact code --- testbeam_analysis/result_analysis.py | 226 +++++++++------------- testbeam_analysis/tools/analysis_utils.py | 50 ++--- 2 files changed, 105 insertions(+), 171 deletions(-) diff --git a/testbeam_analysis/result_analysis.py b/testbeam_analysis/result_analysis.py index b1373a29..680d82fb 100644 --- a/testbeam_analysis/result_analysis.py +++ b/testbeam_analysis/result_analysis.py @@ -1376,6 +1376,12 @@ def normalize_distributions(hists, fit_results): actual_pixel_size = pixel_size[actual_dut] initialize = True + + # create arrays for storing histogrammed intersections for different CS + projections_x_cs_hist = np.empty(shape=(4, n_bins[0])) + projections_y_cs_hist = np.empty(shape=(4, n_bins[1])) + projections_cs_2d_hist = np.empty(shape=(4, n_bins[0], n_bins[1])) + for tracks_chunk, _ in analysis_utils.data_aligned_at_events(node, chunk_size=chunk_size): # read track file in chunks print tracks_chunk.dtype # select only valid track intersection @@ -1384,111 +1390,70 @@ def normalize_distributions(hists, fit_results): # select different CS in order to histogram them and to deduce from the CS ratio the effective CS-1-pitch n_hits = tracks_chunk[~np.isnan(tracks_chunk['n_hits_dut_%d' % actual_dut])]['n_hits_dut_%d' % actual_dut] - # read track intersections wiht actual dut for cluster sizes betwenn 1 and 4 - intersection_x, intersection_y, intersection_z = tracks_chunk[:]['offset_0'], tracks_chunk[:]['offset_1'], tracks_chunk[:]['offset_2'] - intersection_x_cs_1 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] - intersection_y_cs_1 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] - intersection_z_cs_1 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 1)] - intersection_x_cs_2 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] - intersection_y_cs_2 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] - intersection_z_cs_2 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 2)] - intersection_x_cs_3 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] - intersection_y_cs_3 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] - intersection_z_cs_3 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 3)] - intersection_x_cs_4 = intersection_x[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] - intersection_y_cs_4 = intersection_y[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] - intersection_z_cs_4 = intersection_z[np.where(tracks_chunk[:]['n_hits_dut_%i' % actual_dut] == 4)] + intersections = np.column_stack((tracks_chunk[:]['offset_0'], tracks_chunk[:]['offset_1'], tracks_chunk[:]['offset_2'])) + + # arrays for intersections for different CS + intersections_cs_1 = np.zeros(shape=(3, len(n_hits[n_hits == 1]))) + intersections_cs_2 = np.zeros(shape=(3, len(n_hits[n_hits == 2]))) + intersections_cs_3 = np.zeros(shape=(3, len(n_hits[n_hits == 3]))) + intersections_cs_4 = np.zeros(shape=(3, len(n_hits[n_hits == 4]))) + + # read track intersections with actual dut for cluster sizes between 1 and 4 + for dim in range(3): + intersections_cs_1[dim, :] = intersections[:, dim][tracks_chunk['n_hits_dut_%i' % actual_dut] == 1] + intersections_cs_2[dim, :] = intersections[:, dim][tracks_chunk['n_hits_dut_%i' % actual_dut] == 2] + intersections_cs_3[dim, :] = intersections[:, dim][tracks_chunk['n_hits_dut_%i' % actual_dut] == 3] + intersections_cs_4[dim, :] = intersections[:, dim][tracks_chunk['n_hits_dut_%i' % actual_dut] == 4] + + # stack intersections of all CS together + intersections_cs = [intersections_cs_1, intersections_cs_2, intersections_cs_3, intersections_cs_4] + intersections_cs_local = [np.zeros_like(intersections_cs_1), np.zeros_like(intersections_cs_2), + np.zeros_like(intersections_cs_3), np.zeros_like(intersections_cs_4)] # transoform to local coordinate system - if use_prealignment: - intersection_x_cs_1_local, intersection_y_cs_1_local, intersection_z_cs_1_local = geometry_utils.apply_alignment(intersection_x_cs_1, intersection_y_cs_1, intersection_z_cs_1, - dut_index=actual_dut, - prealignment=prealignment, - inverse=True) - - intersection_x_cs_2_local, intersection_y_cs_2_local, intersection_z_cs_2_local = geometry_utils.apply_alignment(intersection_x_cs_2, intersection_y_cs_2, intersection_z_cs_2, - dut_index=actual_dut, - prealignment=prealignment, - inverse=True) - - intersection_x_cs_3_local, intersection_y_cs_3_local, intersection_z_cs_3_local = geometry_utils.apply_alignment(intersection_x_cs_3, intersection_y_cs_3, intersection_z_cs_3, - dut_index=actual_dut, - prealignment=prealignment, - inverse=True) - - intersection_x_cs_4_local, intersection_y_cs_4_local, intersection_z_cs_4_local = geometry_utils.apply_alignment(intersection_x_cs_4, intersection_y_cs_4, intersection_z_cs_4, - dut_index=actual_dut, - prealignment=prealignment, - inverse=True) + for cs in range(4): + if use_prealignment: + intersections_cs_local[cs][0], intersections_cs_local[cs][1], intersections_cs_local[cs][2] = geometry_utils.apply_alignment(intersections_cs[cs][0], intersections_cs[cs][1], intersections_cs[cs][2], + dut_index=actual_dut, + prealignment=prealignment, + inverse=True) + else: + intersections_cs_local[cs][0], intersections_cs_local[cs][1], intersections_cs_local[cs][2] = geometry_utils.apply_alignment(intersections_cs[cs][0], intersections_cs[cs][1], intersections_cs[cs][2], + dut_index=actual_dut, + alignment=alignment, + inverse=True) - else: - intersection_x_cs_1_local, intersection_y_cs_1_local, intersection_z_cs_1_local = geometry_utils.apply_alignment(intersection_x_cs_1, intersection_y_cs_1, intersection_z_cs_1, - dut_index=actual_dut, - alignment=alignment, - inverse=True) - - intersection_x_cs_2_local, intersection_y_cs_2_local, intersection_z_cs_2_local = geometry_utils.apply_alignment(intersection_x_cs_2, intersection_y_cs_2, intersection_z_cs_2, - dut_index=actual_dut, - alignment=alignment, - inverse=True) - - intersection_x_cs_3_local, intersection_y_cs_3_local, intersection_z_cs_3_local = geometry_utils.apply_alignment(intersection_x_cs_3, intersection_y_cs_3, intersection_z_cs_3, - dut_index=actual_dut, - alignment=alignment, - inverse=True) - - intersection_x_cs_4_local, intersection_y_cs_4_local, intersection_z_cs_4_local = geometry_utils.apply_alignment(intersection_x_cs_4, intersection_y_cs_4, intersection_z_cs_4, - dut_index=actual_dut, - alignment=alignment, - inverse=True) - - if not np.allclose(intersection_z_cs_4_local[np.isfinite(intersection_z_cs_4_local)], 0.0): - print intersection_z_cs_4_local - raise RuntimeError("Transformation into local coordinate system gives z != 0") + if not np.allclose(intersections_cs_local[cs][2][np.isfinite(intersections_cs_local[cs][2])], 0.0): + raise RuntimeError("Transformation into local coordinate system gives z != 0") # project track intersections onto one pixel in order to increase statistics - projection_x_1 = np.mod(intersection_x_cs_1_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_1_local))) - projection_y_1 = np.mod(intersection_y_cs_1_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_1_local))) - projection_x_2 = np.mod(intersection_x_cs_2_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_2_local))) - projection_y_2 = np.mod(intersection_y_cs_2_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_2_local))) - projection_x_3 = np.mod(intersection_x_cs_3_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_3_local))) - projection_y_3 = np.mod(intersection_y_cs_3_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_3_local))) - projection_x_4 = np.mod(intersection_x_cs_4_local, np.array([actual_pixel_size[0]] * len(intersection_x_cs_4_local))) - projection_y_4 = np.mod(intersection_y_cs_4_local, np.array([actual_pixel_size[1]] * len(intersection_y_cs_4_local))) - - # histogram intersections and CS (for calculation of effective pitch) + projections_cs = [np.zeros_like(intersections_cs_1), np.zeros_like(intersections_cs_2), + np.zeros_like(intersections_cs_3), np.zeros_like(intersections_cs_4)] + for cs in range(4): + for dim in range(2): + projections_cs[cs][dim] = np.mod(intersections_cs_local[cs][dim], + np.array([actual_pixel_size[dim]] * len(intersections_cs_local[cs][dim]))) + + # histogram intersections and create cluster size histogram (for calculation of effective pitch) if initialize: - projection_x_1_hist, edges_x = np.histogram(projection_x_1, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) - projection_y_1_hist, edges_y = np.histogram(projection_y_1, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) - projection_x_2_hist, _ = np.histogram(projection_x_2, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) - projection_y_2_hist, _ = np.histogram(projection_y_2, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) - projection_x_3_hist, _ = np.histogram(projection_x_3, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) - projection_y_3_hist, _ = np.histogram(projection_y_3, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) - projection_x_4_hist, _ = np.histogram(projection_x_4, bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) - projection_y_4_hist, _ = np.histogram(projection_y_4, bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) - - projection_1_2d_hist, _, _ = np.histogram2d(x=projection_x_1, y=projection_y_1, bins=[edges_x, edges_y]) - projection_2_2d_hist, _, _ = np.histogram2d(x=projection_x_2, y=projection_y_2, bins=[edges_x, edges_y]) - projection_3_2d_hist, _, _ = np.histogram2d(x=projection_x_3, y=projection_y_3, bins=[edges_x, edges_y]) - projection_4_2d_hist, _, _ = np.histogram2d(x=projection_x_4, y=projection_y_4, bins=[edges_x, edges_y]) + for cs in range(4): + if cs == 0: + projections_x_cs_hist[cs], edges_x = np.histogram(projections_cs[cs][0], bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projections_y_cs_hist[cs], edges_y = np.histogram(projections_cs[cs][1], bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + else: + projections_x_cs_hist[cs], _ = np.histogram(projections_cs[cs][0], bins=n_bins[0], range=[0.0, actual_pixel_size[0]]) + projections_y_cs_hist[cs], _ = np.histogram(projections_cs[cs][1], bins=n_bins[1], range=[0.0, actual_pixel_size[1]]) + projections_cs_2d_hist[cs], _, _ = np.histogram2d(x=projections_cs[cs][0], y=projections_cs[cs][1], bins=[edges_x, edges_y]) cs_edges = np.arange(1.0 - 0.5, 1000 + 0.5, 1) cs_hist, cs_edges = np.histogram(n_hits, bins=cs_edges, density=False) initialize = False else: # if already read first chunk, add histograms up - projection_x_1_hist += np.histogram(projection_x_1, bins=edges_x)[0] - projection_y_1_hist += np.histogram(projection_y_1, bins=edges_y)[0] - projection_x_2_hist += np.histogram(projection_x_2, bins=edges_x)[0] - projection_y_2_hist += np.histogram(projection_y_2, bins=edges_y)[0] - projection_x_3_hist += np.histogram(projection_x_3, bins=edges_x)[0] - projection_y_3_hist += np.histogram(projection_y_3, bins=edges_y)[0] - projection_x_4_hist += np.histogram(projection_x_4, bins=edges_x)[0] - projection_y_4_hist += np.histogram(projection_y_4, bins=edges_y)[0] - - projection_1_2d_hist += np.histogram2d(x=projection_x_1, y=projection_y_1, bins=[edges_x, edges_y])[0] - projection_2_2d_hist += np.histogram2d(x=projection_x_2, y=projection_y_2, bins=[edges_x, edges_y])[0] - projection_3_2d_hist += np.histogram2d(x=projection_x_3, y=projection_y_3, bins=[edges_x, edges_y])[0] - projection_4_2d_hist += np.histogram2d(x=projection_x_4, y=projection_y_4, bins=[edges_x, edges_y])[0] + for cs in range(4): + projections_x_cs_hist[cs, :] += np.histogram(projections_cs[cs][0], bins=edges_x)[0] + projections_y_cs_hist[cs, :] += np.histogram(projections_cs[cs][1], bins=edges_y)[0] + + projections_cs_2d_hist[cs] += np.histogram2d(x=projections_cs[cs][0], y=projections_cs[cs][1], bins=[edges_x, edges_y])[0] cs_hist += np.histogram(n_hits, bins=cs_edges, density=False)[0] @@ -1498,34 +1463,14 @@ def normalize_distributions(hists, fit_results): if plot: # plot histograms for cluster sizes between 1 and 4 - plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_1_hist, - y_intersections_hist=projection_y_1_hist, - intersections_2d_hist=projection_1_2d_hist, - pixel_pitch=actual_pixel_size, - output_pdf=output_pdf, - bins=[edges_x, edges_y], - plot_title=('In-pixel Cluster Size 1 Hit Distribution for DUT%d' % actual_dut)) - plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_2_hist, - y_intersections_hist=projection_y_2_hist, - intersections_2d_hist=projection_2_2d_hist, - pixel_pitch=actual_pixel_size, - output_pdf=output_pdf, - bins=[edges_x, edges_y], - plot_title=('In-pixel Cluster Size 2 Hit Distribution for DUT%d' % actual_dut)) - plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_3_hist, - y_intersections_hist=projection_y_3_hist, - intersections_2d_hist=projection_3_2d_hist, - pixel_pitch=actual_pixel_size, - output_pdf=output_pdf, - bins=[edges_x, edges_y], - plot_title=('In-pixel Cluster Size 3 Hit Distribution for DUT%d' % actual_dut)) - plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projection_x_4_hist, - y_intersections_hist=projection_y_4_hist, - intersections_2d_hist=projection_4_2d_hist, - pixel_pitch=actual_pixel_size, - output_pdf=output_pdf, - bins=[edges_x, edges_y], - plot_title=('In-pixel Cluster Size 4 Hit Distribution for DUT%d' % actual_dut)) + for cs in range(4): + plot_utils.plot_in_pixel_hit_hist(x_intersections_hist=projections_x_cs_hist[cs], + y_intersections_hist=projections_y_cs_hist[cs], + intersections_2d_hist=projections_cs_2d_hist[cs], + pixel_pitch=actual_pixel_size, + output_pdf=output_pdf, + bins=[edges_x, edges_y], + plot_title=('In-pixel Cluster Size %i Hit Distribution for DUT%d' % (cs + 1, actual_dut))) # take slice out of middle slice_width = 5 # set withd of slice to 5 um @@ -1537,28 +1482,41 @@ def normalize_distributions(hists, fit_results): stop_x = int(slice_stop_x / (actual_pixel_size[0] / n_bins[0])) start_y = int(slice_start_y / (actual_pixel_size[0] / n_bins[1])) stop_y = int(slice_stop_y / (actual_pixel_size[0] / n_bins[1])) - sliced_hist_y_1 = np.sum(projection_1_2d_hist[start_x:stop_x, ], axis=0) - sliced_hist_x_1 = np.sum(projection_1_2d_hist[:, start_y:stop_y], axis=1) - sliced_hist_y_2 = np.sum(projection_2_2d_hist[start_x:stop_x, ], axis=0) - sliced_hist_x_2 = np.sum(projection_2_2d_hist[:, start_y:stop_y], axis=1) + + # calculate sliced hists + sliced_cs_hists_x = np.empty(shape=(2, n_bins[0])) + sliced_cs_hists_y = np.empty(shape=(2, n_bins[1])) + for cs in range(2): + sliced_cs_hists_y[cs] = np.sum(projections_cs_2d_hist[cs][start_x:stop_x, ], axis=0) + sliced_cs_hists_x[cs] = np.sum(projections_cs_2d_hist[cs][:, start_y:stop_y], axis=1) + + sliced_hist_cs_all = [sliced_cs_hists_x[0], sliced_cs_hists_x[1], sliced_cs_hists_y[0], sliced_cs_hists_y[1]] # fit projections onto x and y-axis of CS 1 and CS 2 distributions - fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced = analysis_utils.fit_in_pixel_hit_hist(cs1_hist_projected=[sliced_hist_x_1, sliced_hist_y_1], - cs2_hist_projected=[sliced_hist_x_2, sliced_hist_y_2], - bin_center_x=(edges_x[1:] + edges_x[:-1]) / 2.0, - bin_center_y=(edges_y[1:] + edges_y[:-1]) / 2.0) + fit_params_cs_sliced = np.zeros(shape=(2, 2, 4)) + for cs in range(2): + fit_params_cs_sliced[cs][0] = analysis_utils.fit_in_pixel_hit_hist(hist=sliced_cs_hists_x[cs], + edges=(edges_x[1:] + edges_x[:-1]) / 2.0) + fit_params_cs_sliced[cs][1] = analysis_utils.fit_in_pixel_hit_hist(hist=sliced_cs_hists_y[cs], + edges=(edges_y[1:] + edges_y[:-1]) / 2.0) + + fit_params_all = [fit_params_cs_sliced[0][0], fit_params_cs_sliced[0][1], fit_params_cs_sliced[1][0], fit_params_cs_sliced[1][1]] + + if np.any(np.isnan(fit_params_all)): + logging.warning("Some fits could not performed.") + continue # deduce normalization (shift) from fit, norm such that background of CS 1 distribution is the same as background of CS 2 distribution - x_intersections_hist_cs_1_normed, x_intersections_hist_cs_2_normed, y_intersections_hist_cs_1_normed, y_intersections_hist_cs_2_normed = normalize_distributions(hists=[sliced_hist_x_1, sliced_hist_x_2, sliced_hist_y_1, sliced_hist_y_2], - fit_results=[fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced]) + x_intersections_hist_cs_1_normed, x_intersections_hist_cs_2_normed, y_intersections_hist_cs_1_normed, y_intersections_hist_cs_2_normed = normalize_distributions(hists=sliced_hist_cs_all, fit_results=fit_params_all) + if plot: # plot effective CS-1 pitch from ratio of CS distribution plot_utils.plot_in_pixel_hit_hist_with_eff_pitch(x_intersections_hist_cs_1=x_intersections_hist_cs_1_normed, x_intersections_hist_cs_2=x_intersections_hist_cs_2_normed, y_intersections_hist_cs_1=y_intersections_hist_cs_1_normed, y_intersections_hist_cs_2=y_intersections_hist_cs_2_normed, - intersections_2d_hist=projection_1_2d_hist - projection_2_2d_hist, - fit_results=[fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced], + intersections_2d_hist=projections_cs_2d_hist[0] - projections_cs_2d_hist[1], + fit_results=fit_params_all, pixel_pitch=actual_pixel_size, bins=[edges_x, edges_y], effective_pitch=effective_pitch, diff --git a/testbeam_analysis/tools/analysis_utils.py b/testbeam_analysis/tools/analysis_utils.py index c02284b6..27beb3da 100644 --- a/testbeam_analysis/tools/analysis_utils.py +++ b/testbeam_analysis/tools/analysis_utils.py @@ -924,43 +924,19 @@ def fit_residuals_vs_position(hist, xedges, yedges, xlabel="", ylabel="", title= return fit, cov -def fit_in_pixel_hit_hist(cs1_hist_projected, cs2_hist_projected, bin_center_x, bin_center_y): - # fit distributions with gauss curve - fit_cs_1_x_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, - bin_center_x, - cs1_hist_projected[0], - p0=(np.max(cs1_hist_projected[0]) - np.min(cs1_hist_projected[0]), - np.median(bin_center_x), - np.std(np.repeat((bin_center_x).astype(np.int64), (cs1_hist_projected[0]).astype(np.int64))), - np.min(cs1_hist_projected[0])), - sigma=np.sqrt(cs1_hist_projected[0])) - fit_cs_1_y_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, - bin_center_y, - cs1_hist_projected[1], - p0=(np.max(cs1_hist_projected[1]) - np.min(cs1_hist_projected[1]), - np.median(bin_center_y), - np.std(np.repeat((bin_center_y).astype(np.int64), (cs1_hist_projected[1]).astype(np.int64))), - np.min(cs1_hist_projected[1])), - sigma=np.sqrt(cs1_hist_projected[1])) - - fit_cs_2_x_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, - bin_center_x, - cs2_hist_projected[0], - p0=(np.max(cs2_hist_projected[0]) - np.min(cs2_hist_projected[0]), - np.median(bin_center_x), - np.std(np.repeat((bin_center_x).astype(np.int64), (cs2_hist_projected[0]).astype(np.int64))), - np.min(cs2_hist_projected[0])), - sigma=np.sqrt(cs2_hist_projected[0])) - fit_cs_2_y_sliced, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, - bin_center_y, - cs2_hist_projected[1], - p0=(np.max(cs2_hist_projected[1]) - np.min(cs2_hist_projected[1]), - np.median(bin_center_y), - np.std(np.repeat((bin_center_y).astype(np.int64), (cs2_hist_projected[1]).astype(np.int64))), - np.min(cs2_hist_projected[1])), - sigma=np.sqrt(cs2_hist_projected[1])) - - return fit_cs_1_x_sliced, fit_cs_1_y_sliced, fit_cs_2_x_sliced, fit_cs_2_y_sliced +def fit_in_pixel_hit_hist(hist, edges): + try: + fit, _ = curve_fit(testbeam_analysis.tools.analysis_utils.gauss_offset, + edges, hist, + p0=(np.max(hist) - np.min(hist), + np.median(edges), + np.std(np.repeat((edges).astype(np.int64), (hist).astype(np.int64))), + np.min(hist)), + sigma=np.sqrt(hist)) + except RuntimeError: + fit = [np.NaN, np.NaN, np.NaN, np.NaN] + + return fit def hough_transform(img, theta_res=1.0, rho_res=1.0, return_edges=False): From a434a1e9f301948241bb0caa823b0bf9c841e1cb Mon Sep 17 00:00:00 2001 From: YannickDieter Date: Thu, 7 Dec 2017 12:25:44 +0100 Subject: [PATCH 3/4] ENH: calculate effective pitch in both dimensions --- testbeam_analysis/result_analysis.py | 7 ++++--- testbeam_analysis/tools/plot_utils.py | 8 ++++---- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/testbeam_analysis/result_analysis.py b/testbeam_analysis/result_analysis.py index 680d82fb..e3fd41c4 100644 --- a/testbeam_analysis/result_analysis.py +++ b/testbeam_analysis/result_analysis.py @@ -1333,10 +1333,11 @@ def calculate_in_pixel_hit_distribution(input_tracks_file, input_alignment_file, use_prealignment = True if force_prealignment else False def calculate_effective_cs1_pitch(cs_hist, pitch): - '''Calculates the effectice CS-1 pitch p_eff using the CS distribution: p_eff = sqrt(pitch^2 * CS1_percentage) + '''Calculates the effectice CS-1 pitch p_eff using the CS distribution: p_eff = sqrt(pitch^2 * N_CS1 / N_tot) + Where N_CS1 is the number of clusters with CS1 and N_tot the total number of clusters. ''' perc_cluster_size_1 = 100. * ((cs_hist[0].astype(np.float) / np.sum(cs_hist).astype(np.float))) - return np.sqrt(pitch * pitch * perc_cluster_size_1 / 100.) + return pitch[0] * np.sqrt(perc_cluster_size_1 / 100.), pitch[1] * np.sqrt(perc_cluster_size_1 / 100.) def normalize_distributions(hists, fit_results): '''Normalize CS distributions in order to compare them. @@ -1458,7 +1459,7 @@ def normalize_distributions(hists, fit_results): cs_hist += np.histogram(n_hits, bins=cs_edges, density=False)[0] # calculate effective pitch from cluster size ratio - effective_pitch = calculate_effective_cs1_pitch(cs_hist=cs_hist, pitch=actual_pixel_size[0]) + effective_pitch = calculate_effective_cs1_pitch(cs_hist=cs_hist, pitch=[actual_pixel_size[0], actual_pixel_size[1]]) effective_pitches.append(effective_pitch) if plot: diff --git a/testbeam_analysis/tools/plot_utils.py b/testbeam_analysis/tools/plot_utils.py index 41534565..4a6defdf 100644 --- a/testbeam_analysis/tools/plot_utils.py +++ b/testbeam_analysis/tools/plot_utils.py @@ -1508,8 +1508,8 @@ def difference_x_sliced(x): def difference_y_sliced(x): return testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_1_y_sliced) - fit_cs_1_y_sliced[3] + fit_cs_2_y_sliced[0] + fit_cs_2_y_sliced[3] - testbeam_analysis.tools.analysis_utils.gauss_offset(x, *fit_cs_2_y_sliced) - roots_x_sliced = fsolve(difference_x_sliced, [(pixel_pitch[0] - effective_pitch) / 2.0, pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.0)]) - roots_y_sliced = fsolve(difference_y_sliced, [(pixel_pitch[1] - effective_pitch) / 2.0, pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.0)]) + roots_x_sliced = fsolve(difference_x_sliced, [(pixel_pitch[0] - effective_pitch[0]) / 2.0, pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch[0]) / 2.0)]) + roots_y_sliced = fsolve(difference_y_sliced, [(pixel_pitch[1] - effective_pitch[1]) / 2.0, pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch[1]) / 2.0)]) # draw intersections of CS 1 and CS 2 distribution for i in range(len(roots_x_sliced)): @@ -1520,8 +1520,8 @@ def difference_y_sliced(x): [roots_y_sliced[0], roots_y_sliced[1], roots_y_sliced[1], roots_y_sliced[0], roots_y_sliced[0]], '-', color='black', label='Effective CS-1-Pitch from Transition of Cluster Size Hit Distributions') # draw effective cluster size area from cluster size ratio - axHeatmap.plot([(pixel_pitch[0] - effective_pitch) / 2., (pixel_pitch[0] - effective_pitch) / 2., pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.), pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch) / 2.), (pixel_pitch[0] - effective_pitch) / 2.], - [(pixel_pitch[1] - effective_pitch) / 2., pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.), pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch) / 2.), (pixel_pitch[1] - effective_pitch) / 2., (pixel_pitch[1] - effective_pitch) / 2], + axHeatmap.plot([(pixel_pitch[0] - effective_pitch[0]) / 2., (pixel_pitch[0] - effective_pitch[0]) / 2., pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch[0]) / 2.), pixel_pitch[0] - ((pixel_pitch[0] - effective_pitch[0]) / 2.), (pixel_pitch[0] - effective_pitch[0]) / 2.], + [(pixel_pitch[1] - effective_pitch[1]) / 2., pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch[1]) / 2.), pixel_pitch[1] - ((pixel_pitch[1] - effective_pitch[1]) / 2.), (pixel_pitch[1] - effective_pitch[1]) / 2., (pixel_pitch[1] - effective_pitch[1]) / 2], '--', color='black', label='Effective CS-1-Pitch from Cluster Size Ratio of Cluster Size Distribution') # colorbar From 94bdd2d32aff923b018fe518076f44d55e1446eb Mon Sep 17 00:00:00 2001 From: YannickDieter Date: Thu, 7 Dec 2017 12:26:21 +0100 Subject: [PATCH 4/4] MAINT: clean up --- testbeam_analysis/result_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testbeam_analysis/result_analysis.py b/testbeam_analysis/result_analysis.py index e3fd41c4..c3bdd217 100644 --- a/testbeam_analysis/result_analysis.py +++ b/testbeam_analysis/result_analysis.py @@ -1384,7 +1384,6 @@ def normalize_distributions(hists, fit_results): projections_cs_2d_hist = np.empty(shape=(4, n_bins[0], n_bins[1])) for tracks_chunk, _ in analysis_utils.data_aligned_at_events(node, chunk_size=chunk_size): # read track file in chunks - print tracks_chunk.dtype # select only valid track intersection tracks_chunk = tracks_chunk[~np.isnan(tracks_chunk[:]['offset_0'])] @@ -1449,6 +1448,7 @@ def normalize_distributions(hists, fit_results): cs_edges = np.arange(1.0 - 0.5, 1000 + 0.5, 1) cs_hist, cs_edges = np.histogram(n_hits, bins=cs_edges, density=False) initialize = False + else: # if already read first chunk, add histograms up for cs in range(4): projections_x_cs_hist[cs, :] += np.histogram(projections_cs[cs][0], bins=edges_x)[0]