<commit_before>''' All functions creating results (e.g. efficiency, residuals, track density) from fitted tracks are listed here.'''
from __future__ import division

import logging
import re
from collections import Iterable
import os.path

import tables as tb
import numpy as np
from matplotlib.backends.backend_pdf import PdfPages
from scipy.stats import binned_statistic_2d
from scipy.optimize import curve_fit

from testbeam_analysis.tools import plot_utils
from testbeam_analysis.tools import geometry_utils
from testbeam_analysis.tools import analysis_utils


def calculate_residuals(input_tracks_file, input_alignment_file, n_pixels, pixel_size, output_residuals_file=None, dut_names=None, use_duts=None, max_chi2=None, nbins_per_pixel=None, npixels_per_bin=None, force_prealignment=False, plot=True, chunk_size=1000000):
    '''Takes the tracks and calculates residuals for selected DUTs in col, row direction.

    Parameters
    ----------
    input_tracks_file : string
        Filename of the input tracks file.
    input_alignment_file : string
        Filename of the input aligment file.
    n_pixels : iterable of tuples
        One tuple per DUT describing the number of pixels in column, row direction
        e.g. for 2 DUTs: n_pixels = [(80, 336), (80, 336)]
    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)]
    output_residuals_file : string
        Filename of the output residuals file. If None, the filename will be derived from the input hits file.
    dut_names : iterable
        Name of the DUTs. If None, DUT numbers will be used.
    use_duts : iterable
        The duts to calculate residuals for. If None all duts in the input_tracks_file are used
    max_chi2 : uint, iterable
        Use only not heavily scattered tracks to increase track pointing resolution (cut on chi2).
        Cut can be a number and is used then for all DUTS or a list with a chi 2 cut for each DUT.
        If None, no cut is applied.
    nbins_per_pixel : int
        Number of bins per pixel along the residual axis. Number is a positive integer or None to automatically set the binning.
    npixels_per_bin : int
        Number of pixels per bin along the position axis. Number is a positive integer or None to automatically set the binning.
    force_prealignment : boolean
        Take the prealignment, although if a coarse alignment is availale.
    plot : bool
        If True, create additional output plots.
    chunk_size : int
        Chunk size of the data when reading from file.
    '''
    logging.info('=== Calculating residuals ===')

    use_prealignment = True if force_prealignment else False

    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]

    if output_residuals_file is None:
        output_residuals_file = os.path.splitext(input_tracks_file)[0] + '_residuals.h5'

    if plot is True:
        output_pdf = PdfPages(os.path.splitext(output_residuals_file)[0] + '.pdf')
    else:
        output_pdf = None

    if not isinstance(max_chi2, Iterable):
        max_chi2 = [max_chi2] * n_duts

    with tb.open_file(input_tracks_file, mode='r') as in_file_h5:
        with tb.open_file(output_residuals_file, mode='w') as out_file_h5:
            for node in in_file_h5.root:
                actual_dut = int(re.findall(r'\d+', node.name)[-1])
                if use_duts and actual_dut not in use_duts:
                    continue
                logging.debug('Calculate residuals for DUT%d', actual_dut)

                initialize = True  # initialize the histograms
                for tracks_chunk, _ in analysis_utils.data_aligned_at_events(node, chunk_size=chunk_size):
                    # select good hits and tracks
                    selection = np.logical_and(~np.isnan(tracks_chunk['x_dut_%d' % actual_dut]), ~np.isnan(tracks_chunk['track_chi2']))
                    tracks_chunk = tracks_chunk[selection]  # Take only tracks where actual dut has a hit, otherwise residual wrong
                    if max_chi2[actual_dut] is not None:
                        tracks_chunk = tracks_chunk[tracks_chunk['track_chi2'] <= max_chi2[actual_dut]]

                    # Coordinates in global coordinate system (x, y, z)
                    hit_x, hit_y, hit_z = tracks_chunk['x_dut_%d' % actual_dut], tracks_chunk['y_dut_%d' % actual_dut], tracks_chunk['z_dut_%d' % actual_dut]
                    intersection_x, intersection_y, intersection_z = tracks_chunk['offset_0'], tracks_chunk['offset_1'], tracks_chunk['offset_2']

                    # Transform to local coordinate system
                    if use_prealignment:
                        hit_x_local, hit_y_local, hit_z_local = geometry_utils.apply_alignment(hit_x, hit_y, hit_z,
                                                                                               dut_index=actual_dut,
                                                                                               prealignment=prealignment,
                                                                                               inverse=True)
                        intersection_x_local, intersection_y_local, intersection_z_local = geometry_utils.apply_alignment(intersection_x, intersection_y, intersection_z,
                                                                                                                          dut_index=actual_dut,
                                                                                                                          prealignment=prealignment,
                                                                                                                          inverse=True)
                    else:  # Apply transformation from fine alignment information
                        hit_x_local, hit_y_local, hit_z_local = geometry_utils.apply_alignment(hit_x, hit_y, hit_z,
                                                                                               dut_index=actual_dut,
                                                                                               alignment=alignment,
                                                                                               inverse=True)
                        intersection_x_local, intersection_y_local, intersection_z_local = geometry_utils.apply_alignment(intersection_x, intersection_y, intersection_z,
                                                                                                                          dut_index=actual_dut,
                                                                                                                          alignment=alignment,
                                                                                                                          inverse=True)

                    if not np.allclose(hit_z_local, 0.0) or not np.allclose(intersection_z_local, 0.0):
                        logging.error('Hit z position = %s and z intersection %s', str(hit_z_local[:3]), str(intersection_z_local[:3]))
                        raise RuntimeError('The transformation to the local coordinate system did not give all z = 0. Wrong alignment used?')

                    difference = np.column_stack((hit_x, hit_y, hit_z)) - np.column_stack((intersection_x, intersection_y, intersection_z))
                    difference_local = np.column_stack((hit_x_local, hit_y_local, hit_z_local)) - np.column_stack((intersection_x_local, intersection_y_local, intersection_z_local))

                    # Histogram residuals in different ways
                    if initialize:  # Only true for the first iteration, calculate the binning for the histograms
                        initialize = False

                        # detect peaks and calculate width to estimate the size of the histograms
                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference[:, 0]), np.max(difference[:, 0])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][0] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][0] / nbins_per_pixel), pixel_size[actual_dut][0] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference[:, 0], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_x, fwhm_x, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_x, fwhm_x, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference[:, 1]), np.max(difference[:, 1])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][1] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][1] / nbins_per_pixel), pixel_size[actual_dut][1] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference[:, 1], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_y, fwhm_y, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_y, fwhm_y, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference_local[:, 0]), np.max(difference_local[:, 0])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][0] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][0] / nbins_per_pixel), pixel_size[actual_dut][0] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference_local[:, 0], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_col, fwhm_col, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_col, fwhm_col, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference_local[:, 1]), np.max(difference_local[:, 1])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][1] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][1] / nbins_per_pixel), pixel_size[actual_dut][1] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference_local[:, 1], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_row, fwhm_row, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_row, fwhm_row, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        # calculate the binning of the histograms, the minimum size is given by plot_npixels, otherwise FWHM is taken into account
                        plot_npixels = 6.0
                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][0], pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][0]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][0])
                            x_range = (center_x - 0.5 * width, center_x + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0])
                            x_range = (center_x - width, center_x + width)
                        hist_residual_x_hist, hist_residual_x_xedges = np.histogram(difference[:, 0], range=x_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_x), np.max(intersection_x)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][1], npixels_per_bin * pixel_size[actual_dut][1])
                        else:
                            nbins = "auto"
                        _, hist_residual_x_yedges = np.histogram(intersection_x, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][1], pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][1]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][1]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][1])
                            y_range = (center_y - 0.5 * width, center_y + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_y / pixel_size[actual_dut][1])
                            y_range = (center_y - width, center_y + width)
                        hist_residual_y_hist, hist_residual_y_yedges = np.histogram(difference[:, 1], range=y_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_y), np.max(intersection_y)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][0], npixels_per_bin * pixel_size[actual_dut][0])
                        else:
                            nbins = "auto"
                        _, hist_residual_y_xedges = np.histogram(intersection_y, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][0], pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][0]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][0])
                            col_range = (center_col - 0.5 * width, center_col + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_col / pixel_size[actual_dut][0])
                            col_range = (center_col - width, center_col + width)
                        hist_residual_col_hist, hist_residual_col_xedges = np.histogram(difference_local[:, 0], range=col_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_x_local), np.max(intersection_x_local)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][1], npixels_per_bin * pixel_size[actual_dut][1])
                        else:
                            nbins = "auto"
                        _, hist_residual_col_yedges = np.histogram(intersection_x_local, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][1], pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][1]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][1]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][1])
                            row_range = (center_row - 0.5 * width, center_row + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_row / pixel_size[actual_dut][1])
                            row_range = (center_row - width, center_row + width)
                        hist_residual_row_hist, hist_residual_row_yedges = np.histogram(difference_local[:, 1], range=row_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_y_local), np.max(intersection_y_local)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][0], npixels_per_bin * pixel_size[actual_dut][0])
                        else:
                            nbins = "auto"
                        _, hist_residual_row_xedges = np.histogram(intersection_y_local, bins=nbins)

                        # global x residual against x position
                        hist_x_residual_x_hist, hist_x_residual_x_xedges, hist_x_residual_x_yedges = np.histogram2d(
                            intersection_x,
                            difference[:, 0],
                            bins=(hist_residual_x_yedges, hist_residual_x_xedges))

                        # global y residual against y position
                        hist_y_residual_y_hist, hist_y_residual_y_xedges, hist_y_residual_y_yedges = np.histogram2d(
                            intersection_y,
                            difference[:, 1],
                            bins=(hist_residual_y_xedges, hist_residual_y_yedges))

                        # global y residual against x position
                        hist_x_residual_y_hist, hist_x_residual_y_xedges, hist_x_residual_y_yedges = np.histogram2d(
                            intersection_x,
                            difference[:, 1],
                            bins=(hist_residual_x_yedges, hist_residual_y_yedges))

                        # global x residual against y position
                        hist_y_residual_x_hist, hist_y_residual_x_xedges, hist_y_residual_x_yedges = np.histogram2d(
                            intersection_y,
                            difference[:, 0],
                            bins=(hist_residual_y_xedges, hist_residual_x_xedges))

                        # local column residual against column position
                        hist_col_residual_col_hist, hist_col_residual_col_xedges, hist_col_residual_col_yedges = np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 0],
                            bins=(hist_residual_col_yedges, hist_residual_col_xedges))

                        # local row residual against row position
                        hist_row_residual_row_hist, hist_row_residual_row_xedges, hist_row_residual_row_yedges = np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 1],
                            bins=(hist_residual_row_xedges, hist_residual_row_yedges))

                        # local row residual against column position
                        hist_col_residual_row_hist, hist_col_residual_row_xedges, hist_col_residual_row_yedges = np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 1],
                            bins=(hist_residual_col_yedges, hist_residual_row_yedges))

                        # local column residual against row position
                        hist_row_residual_col_hist, hist_row_residual_col_xedges, hist_row_residual_col_yedges = np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 0],
                            bins=(hist_residual_row_xedges, hist_residual_col_xedges))

                    else:  # adding data to existing histograms
                        hist_residual_x_hist += np.histogram(difference[:, 0], bins=hist_residual_x_xedges)[0]
                        hist_residual_y_hist += np.histogram(difference[:, 1], bins=hist_residual_y_yedges)[0]
                        hist_residual_col_hist += np.histogram(difference_local[:, 0], bins=hist_residual_col_xedges)[0]
                        hist_residual_row_hist += np.histogram(difference_local[:, 1], bins=hist_residual_row_yedges)[0]

                        # global x residual against x position
                        hist_x_residual_x_hist += np.histogram2d(
                            intersection_x,
                            difference[:, 0],
                            bins=(hist_x_residual_x_xedges, hist_x_residual_x_yedges))[0]

                        # global y residual against y position
                        hist_y_residual_y_hist += np.histogram2d(
                            intersection_y,
                            difference[:, 1],
                            bins=(hist_y_residual_y_xedges, hist_y_residual_y_yedges))[0]

                        # global y residual against x position
                        hist_x_residual_y_hist += np.histogram2d(
                            intersection_x,
                            difference[:, 1],
                            bins=(hist_x_residual_y_xedges, hist_x_residual_y_yedges))[0]

                        # global x residual against y position
                        hist_y_residual_x_hist += np.histogram2d(
                            intersection_y,
                            difference[:, 0],
                            bins=(hist_y_residual_x_xedges, hist_y_residual_x_yedges))[0]

                        # local column residual against column position
                        hist_col_residual_col_hist += np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 0],
                            bins=(hist_col_residual_col_xedges, hist_col_residual_col_yedges))[0]

                        # local row residual against row position
                        hist_row_residual_row_hist += np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 1],
                            bins=(hist_row_residual_row_xedges, hist_row_residual_row_yedges))[0]

                        # local row residual against column position
                        hist_col_residual_row_hist += np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 1],
                            bins=(hist_col_residual_row_xedges, hist_col_residual_row_yedges))[0]

                        # local column residual against row position
                        hist_row_residual_col_hist += np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 0],
                            bins=(hist_row_residual_col_xedges, hist_row_residual_col_yedges))[0]

                logging.debug('Storing residual histograms...')

                dut_name = dut_names[actual_dut] if dut_names else ("DUT" + str(actual_dut))
                # Global residuals
                fit_residual_x, cov_residual_x = analysis_utils.fit_residuals(
                    hist=hist_residual_x_hist,
                    edges=hist_residual_x_xedges,
                    label='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                      name='ResidualsX_DUT%d' % (actual_dut),
                                                      title='Residual distribution in x direction for %s' % (dut_name),
                                                      atom=tb.Atom.from_dtype(hist_residual_x_hist.dtype),
                                                      shape=hist_residual_x_hist.shape,
                                                      filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_x.attrs.xedges = hist_residual_x_xedges
                out_res_x.attrs.fit_coeff = fit_residual_x
                out_res_x.attrs.fit_cov = cov_residual_x
                out_res_x[:] = hist_residual_x_hist

                fit_residual_y, cov_residual_y = analysis_utils.fit_residuals(
                    hist=hist_residual_y_hist,
                    edges=hist_residual_y_yedges,
                    label='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                      name='ResidualsY_DUT%d' % (actual_dut),
                                                      title='Residual distribution in y direction for %s' % (dut_name),
                                                      atom=tb.Atom.from_dtype(hist_residual_y_hist.dtype),
                                                      shape=hist_residual_y_hist.shape,
                                                      filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_y.attrs.yedges = hist_residual_y_yedges
                out_res_y.attrs.fit_coeff = fit_residual_y
                out_res_y.attrs.fit_cov = cov_residual_y
                out_res_y[:] = hist_residual_y_hist

                fit_x_residual_x, cov_x_residual_x = analysis_utils.fit_residuals_vs_position(
                    hist=hist_x_residual_x_hist,
                    xedges=hist_x_residual_x_xedges,
                    yedges=hist_x_residual_x_yedges,
                    xlabel='X position [um]',
                    ylabel='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )

                out_x_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                        name='XResidualsX_DUT%d' % (actual_dut),
                                                        title='Residual distribution in x direction as a function of the x position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_x_residual_x_hist.dtype),
                                                        shape=hist_x_residual_x_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_x_res_x.attrs.xedges = hist_x_residual_x_xedges
                out_x_res_x.attrs.yedges = hist_x_residual_x_yedges
                out_x_res_x.attrs.fit_coeff = fit_x_residual_x
                out_x_res_x.attrs.fit_cov = cov_x_residual_x
                out_x_res_x[:] = hist_x_residual_x_hist

                fit_y_residual_y, cov_y_residual_y = analysis_utils.fit_residuals_vs_position(
                    hist=hist_y_residual_y_hist,
                    xedges=hist_y_residual_y_xedges,
                    yedges=hist_y_residual_y_yedges,
                    xlabel='Y position [um]',
                    ylabel='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_y_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                        name='YResidualsY_DUT%d' % (actual_dut),
                                                        title='Residual distribution in y direction as a function of the y position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_y_residual_y_hist.dtype),
                                                        shape=hist_y_residual_y_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_y_res_y.attrs.xedges = hist_y_residual_y_xedges
                out_y_res_y.attrs.yedges = hist_y_residual_y_yedges
                out_y_res_y.attrs.fit_coeff = fit_y_residual_y
                out_y_res_y.attrs.fit_cov = cov_y_residual_y
                out_y_res_y[:] = hist_y_residual_y_hist

                fit_x_residual_y, cov_x_residual_y = analysis_utils.fit_residuals_vs_position(
                    hist=hist_x_residual_y_hist,
                    xedges=hist_x_residual_y_xedges,
                    yedges=hist_x_residual_y_yedges,
                    xlabel='X position [um]',
                    ylabel='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_x_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                        name='XResidualsY_DUT%d' % (actual_dut),
                                                        title='Residual distribution in y direction as a function of the x position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_x_residual_y_hist.dtype),
                                                        shape=hist_x_residual_y_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_x_res_y.attrs.xedges = hist_x_residual_y_xedges
                out_x_res_y.attrs.yedges = hist_x_residual_y_yedges
                out_x_res_y.attrs.fit_coeff = fit_x_residual_y
                out_x_res_y.attrs.fit_cov = cov_x_residual_y
                out_x_res_y[:] = hist_x_residual_y_hist

                fit_y_residual_x, cov_y_residual_x = analysis_utils.fit_residuals_vs_position(
                    hist=hist_y_residual_x_hist,
                    xedges=hist_y_residual_x_xedges,
                    yedges=hist_y_residual_x_yedges,
                    xlabel='Y position [um]',
                    ylabel='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_y_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                        name='YResidualsX_DUT%d' % (actual_dut),
                                                        title='Residual distribution in x direction as a function of the y position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_y_residual_x_hist.dtype),
                                                        shape=hist_y_residual_x_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_y_res_x.attrs.xedges = hist_y_residual_x_xedges
                out_y_res_x.attrs.yedges = hist_y_residual_x_yedges
                out_y_res_x.attrs.fit_coeff = fit_y_residual_x
                out_y_res_x.attrs.fit_cov = cov_y_residual_x
                out_y_res_x[:] = hist_y_residual_x_hist

                # Local residuals
                fit_residual_col, cov_residual_col = analysis_utils.fit_residuals(
                    hist=hist_residual_col_hist,
                    edges=hist_residual_col_xedges,
                    label='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                        name='ResidualsCol_DUT%d' % (actual_dut),
                                                        title='Residual distribution in column direction for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_residual_col_hist.dtype),
                                                        shape=hist_residual_col_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_col.attrs.xedges = hist_residual_col_xedges
                out_res_col.attrs.fit_coeff = fit_residual_col
                out_res_col.attrs.fit_cov = cov_residual_col
                out_res_col[:] = hist_residual_col_hist

                fit_residual_row, cov_residual_row = analysis_utils.fit_residuals(
                    hist=hist_residual_row_hist,
                    edges=hist_residual_row_yedges,
                    label='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                        name='ResidualsRow_DUT%d' % (actual_dut),
                                                        title='Residual distribution in row direction for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_residual_row_hist.dtype),
                                                        shape=hist_residual_row_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_row.attrs.yedges = hist_residual_row_yedges
                out_res_row.attrs.fit_coeff = fit_residual_row
                out_res_row.attrs.fit_cov = cov_residual_row
                out_res_row[:] = hist_residual_row_hist

                fit_col_residual_col, cov_col_residual_col = analysis_utils.fit_residuals_vs_position(
                    hist=hist_col_residual_col_hist,
                    xedges=hist_col_residual_col_xedges,
                    yedges=hist_col_residual_col_yedges,
                    xlabel='Column position [um]',
                    ylabel='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_col_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                            name='ColResidualsCol_DUT%d' % (actual_dut),
                                                            title='Residual distribution in column direction as a function of the column position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_col_residual_col_hist.dtype),
                                                            shape=hist_col_residual_col_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_col_res_col.attrs.xedges = hist_col_residual_col_xedges
                out_col_res_col.attrs.yedges = hist_col_residual_col_yedges
                out_col_res_col.attrs.fit_coeff = fit_col_residual_col
                out_col_res_col.attrs.fit_cov = cov_col_residual_col
                out_col_res_col[:] = hist_col_residual_col_hist

                fit_row_residual_row, cov_row_residual_row = analysis_utils.fit_residuals_vs_position(
                    hist=hist_row_residual_row_hist,
                    xedges=hist_row_residual_row_xedges,
                    yedges=hist_row_residual_row_yedges,
                    xlabel='Row position [um]',
                    ylabel='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_row_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                            name='RowResidualsRow_DUT%d' % (actual_dut),
                                                            title='Residual distribution in row direction as a function of the row position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_row_residual_row_hist.dtype),
                                                            shape=hist_row_residual_row_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_row_res_row.attrs.xedges = hist_row_residual_row_xedges
                out_row_res_row.attrs.yedges = hist_row_residual_row_yedges
                out_row_res_row.attrs.fit_coeff = fit_row_residual_row
                out_row_res_row.attrs.fit_cov = cov_row_residual_row
                out_row_res_row[:] = hist_row_residual_row_hist

                fit_col_residual_row, cov_col_residual_row = analysis_utils.fit_residuals_vs_position(
                    hist=hist_col_residual_row_hist,
                    xedges=hist_col_residual_row_xedges,
                    yedges=hist_col_residual_row_yedges,
                    xlabel='Column position [um]',
                    ylabel='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_col_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                            name='ColResidualsRow_DUT%d' % (actual_dut),
                                                            title='Residual distribution in row direction as a function of the column position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_col_residual_row_hist.dtype),
                                                            shape=hist_col_residual_row_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_col_res_row.attrs.xedges = hist_col_residual_row_xedges
                out_col_res_row.attrs.yedges = hist_col_residual_row_yedges
                out_col_res_row.attrs.fit_coeff = fit_col_residual_row
                out_col_res_row.attrs.fit_cov = cov_col_residual_row
                out_col_res_row[:] = hist_col_residual_row_hist

                fit_row_residual_col, cov_row_residual_col = analysis_utils.fit_residuals_vs_position(
                    hist=hist_row_residual_col_hist,
                    xedges=hist_row_residual_col_xedges,
                    yedges=hist_row_residual_col_yedges,
                    xlabel='Row position [um]',
                    ylabel='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_row_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                            name='RowResidualsCol_DUT%d' % (actual_dut),
                                                            title='Residual distribution in column direction as a function of the row position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_row_residual_col_hist.dtype),
                                                            shape=hist_row_residual_col_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_row_res_col.attrs.xedges = hist_row_residual_col_xedges
                out_row_res_col.attrs.yedges = hist_row_residual_col_yedges
                out_row_res_col.attrs.fit_coeff = fit_row_residual_col
                out_row_res_col.attrs.fit_cov = cov_row_residual_col
                out_row_res_col[:] = hist_row_residual_col_hist

    if output_pdf is not None:
        output_pdf.close()


def calculate_efficiency(input_tracks_file, input_alignment_file, bin_size, sensor_size, output_efficiency_file=None, pixel_size=None, n_pixels=None, minimum_track_density=1, max_distance=500, use_duts=None, max_chi2=None, force_prealignment=False, cut_distance=None, col_range=None, row_range=None, show_inefficient_events=False, plot=True, chunk_size=1000000):
    '''Takes the tracks and calculates the hit efficiency and hit/track hit distance for selected DUTs.

    Parameters
    ----------
    input_tracks_file : string
        Filename of the input tracks file.
    input_alignment_file : string
        Filename of the input alignment file.
    bin_size : iterable
        Sizes of bins (i.e. (virtual) pixel size). Give one tuple (x, y) for every plane or list of tuples for different planes.
    sensor_size : Tuple or list of tuples
        Describes the sensor size for each DUT. If one tuple is given it is (size x, size y)
        If several tuples are given it is [(DUT0 size x, DUT0 size y), (DUT1 size x, DUT1 size y), ...]
    output_efficiency_file : string
        Filename of the output efficiency file. If None, the filename will be derived from the input hits file.
    minimum_track_density : int
        Minimum track density required to consider bin for efficiency calculation.
    use_duts : iterable
        Calculate the efficiency for selected DUTs. If None, all duts are selected.
    max_chi2 : uint
        Only use tracks with a chi2 <= max_chi2.
    force_prealignment : bool
        Take the prealignment, although if a coarse alignment is availale.
    cut_distance : int
        Use only distances (between DUT hit and track hit) smaller than cut_distance.
    max_distance : int
        Defines binnig of distance values.
    col_range, row_range : iterable
        Column / row value to calculate efficiency for (to neglect noisy edge pixels for efficiency calculation).
    plot : bool
        If True, create additional output plots.
    chunk_size : int
        Chunk size of the data when reading from file.
    '''
    logging.info('=== Calculating efficiency ===')

    if output_efficiency_file is None:
        output_efficiency_file = os.path.splitext(input_tracks_file)[0] + '_efficiency.h5'

    if plot is True:
        output_pdf = PdfPages(os.path.splitext(output_efficiency_file)[0] + '.pdf')
    else:
        output_pdf = None

    use_prealignment = True if force_prealignment else False

    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]

    if not isinstance(max_chi2, Iterable):
        max_chi2 = [max_chi2] * n_duts

    efficiencies = []
    pass_tracks = []
    total_tracks = []
    with tb.open_file(input_tracks_file, mode='r') as in_file_h5:
        for index, node in enumerate(in_file_h5.root):
            actual_dut = int(re.findall(r'\d+', node.name)[-1])
            if use_duts and actual_dut not in use_duts:
                continue
            logging.info('Calculate efficiency for DUT%d', actual_dut)<commit_msg>MAINT: fix length of iterable<commit_after>''' All functions creating results (e.g. efficiency, residuals, track density) from fitted tracks are listed here.'''
from __future__ import division

import logging
import re
from collections import Iterable
import os.path

import tables as tb
import numpy as np
from matplotlib.backends.backend_pdf import PdfPages
from scipy.stats import binned_statistic_2d
from scipy.optimize import curve_fit

from testbeam_analysis.tools import plot_utils
from testbeam_analysis.tools import geometry_utils
from testbeam_analysis.tools import analysis_utils


def calculate_residuals(input_tracks_file, input_alignment_file, n_pixels, pixel_size, output_residuals_file=None, dut_names=None, use_duts=None, max_chi2=None, nbins_per_pixel=None, npixels_per_bin=None, force_prealignment=False, plot=True, chunk_size=1000000):
    '''Takes the tracks and calculates residuals for selected DUTs in col, row direction.

    Parameters
    ----------
    input_tracks_file : string
        Filename of the input tracks file.
    input_alignment_file : string
        Filename of the input aligment file.
    n_pixels : iterable of tuples
        One tuple per DUT describing the number of pixels in column, row direction
        e.g. for 2 DUTs: n_pixels = [(80, 336), (80, 336)]
    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)]
    output_residuals_file : string
        Filename of the output residuals file. If None, the filename will be derived from the input hits file.
    dut_names : iterable
        Name of the DUTs. If None, DUT numbers will be used.
    use_duts : iterable
        The duts to calculate residuals for. If None all duts in the input_tracks_file are used
    max_chi2 : uint, iterable
        Use only not heavily scattered tracks to increase track pointing resolution (cut on chi2).
        Cut can be a number and is used then for all DUTS or a list with a chi 2 cut for each DUT.
        If None, no cut is applied.
    nbins_per_pixel : int
        Number of bins per pixel along the residual axis. Number is a positive integer or None to automatically set the binning.
    npixels_per_bin : int
        Number of pixels per bin along the position axis. Number is a positive integer or None to automatically set the binning.
    force_prealignment : boolean
        Take the prealignment, although if a coarse alignment is availale.
    plot : bool
        If True, create additional output plots.
    chunk_size : int
        Chunk size of the data when reading from file.
    '''
    logging.info('=== Calculating residuals ===')

    use_prealignment = True if force_prealignment else False

    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]

    if output_residuals_file is None:
        output_residuals_file = os.path.splitext(input_tracks_file)[0] + '_residuals.h5'

    if plot is True:
        output_pdf = PdfPages(os.path.splitext(output_residuals_file)[0] + '.pdf')
    else:
        output_pdf = None

    if not isinstance(max_chi2, Iterable):
        max_chi2 = [max_chi2] * n_duts

    with tb.open_file(input_tracks_file, mode='r') as in_file_h5:
        with tb.open_file(output_residuals_file, mode='w') as out_file_h5:
            for node in in_file_h5.root:
                actual_dut = int(re.findall(r'\d+', node.name)[-1])
                if use_duts and actual_dut not in use_duts:
                    continue
                logging.debug('Calculate residuals for DUT%d', actual_dut)

                initialize = True  # initialize the histograms
                for tracks_chunk, _ in analysis_utils.data_aligned_at_events(node, chunk_size=chunk_size):
                    # select good hits and tracks
                    selection = np.logical_and(~np.isnan(tracks_chunk['x_dut_%d' % actual_dut]), ~np.isnan(tracks_chunk['track_chi2']))
                    tracks_chunk = tracks_chunk[selection]  # Take only tracks where actual dut has a hit, otherwise residual wrong
                    if max_chi2[actual_dut] is not None:
                        tracks_chunk = tracks_chunk[tracks_chunk['track_chi2'] <= max_chi2[actual_dut]]

                    # Coordinates in global coordinate system (x, y, z)
                    hit_x, hit_y, hit_z = tracks_chunk['x_dut_%d' % actual_dut], tracks_chunk['y_dut_%d' % actual_dut], tracks_chunk['z_dut_%d' % actual_dut]
                    intersection_x, intersection_y, intersection_z = tracks_chunk['offset_0'], tracks_chunk['offset_1'], tracks_chunk['offset_2']

                    # Transform to local coordinate system
                    if use_prealignment:
                        hit_x_local, hit_y_local, hit_z_local = geometry_utils.apply_alignment(hit_x, hit_y, hit_z,
                                                                                               dut_index=actual_dut,
                                                                                               prealignment=prealignment,
                                                                                               inverse=True)
                        intersection_x_local, intersection_y_local, intersection_z_local = geometry_utils.apply_alignment(intersection_x, intersection_y, intersection_z,
                                                                                                                          dut_index=actual_dut,
                                                                                                                          prealignment=prealignment,
                                                                                                                          inverse=True)
                    else:  # Apply transformation from fine alignment information
                        hit_x_local, hit_y_local, hit_z_local = geometry_utils.apply_alignment(hit_x, hit_y, hit_z,
                                                                                               dut_index=actual_dut,
                                                                                               alignment=alignment,
                                                                                               inverse=True)
                        intersection_x_local, intersection_y_local, intersection_z_local = geometry_utils.apply_alignment(intersection_x, intersection_y, intersection_z,
                                                                                                                          dut_index=actual_dut,
                                                                                                                          alignment=alignment,
                                                                                                                          inverse=True)

                    if not np.allclose(hit_z_local, 0.0) or not np.allclose(intersection_z_local, 0.0):
                        logging.error('Hit z position = %s and z intersection %s', str(hit_z_local[:3]), str(intersection_z_local[:3]))
                        raise RuntimeError('The transformation to the local coordinate system did not give all z = 0. Wrong alignment used?')

                    difference = np.column_stack((hit_x, hit_y, hit_z)) - np.column_stack((intersection_x, intersection_y, intersection_z))
                    difference_local = np.column_stack((hit_x_local, hit_y_local, hit_z_local)) - np.column_stack((intersection_x_local, intersection_y_local, intersection_z_local))

                    # Histogram residuals in different ways
                    if initialize:  # Only true for the first iteration, calculate the binning for the histograms
                        initialize = False

                        # detect peaks and calculate width to estimate the size of the histograms
                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference[:, 0]), np.max(difference[:, 0])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][0] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][0] / nbins_per_pixel), pixel_size[actual_dut][0] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference[:, 0], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_x, fwhm_x, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_x, fwhm_x, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference[:, 1]), np.max(difference[:, 1])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][1] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][1] / nbins_per_pixel), pixel_size[actual_dut][1] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference[:, 1], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_y, fwhm_y, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_y, fwhm_y, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference_local[:, 0]), np.max(difference_local[:, 0])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][0] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][0] / nbins_per_pixel), pixel_size[actual_dut][0] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference_local[:, 0], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_col, fwhm_col, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_col, fwhm_col, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        if nbins_per_pixel is not None:
                            min_difference, max_difference = np.min(difference_local[:, 1]), np.max(difference_local[:, 1])
                            nbins = np.arange(min_difference - (pixel_size[actual_dut][1] / nbins_per_pixel), max_difference + 2 * (pixel_size[actual_dut][1] / nbins_per_pixel), pixel_size[actual_dut][1] / nbins_per_pixel)
                        else:
                            nbins = "auto"
                        hist, edges = np.histogram(difference_local[:, 1], bins=nbins)
                        edge_center = (edges[1:] + edges[:-1]) / 2.0
                        try:
                            _, center_row, fwhm_row, _ = analysis_utils.peak_detect(edge_center, hist)
                        except RuntimeError:
                            # do some simple FWHM with numpy array
                            _, center_row, fwhm_row, _ = analysis_utils.simple_peak_detect(edge_center, hist)

                        # calculate the binning of the histograms, the minimum size is given by plot_npixels, otherwise FWHM is taken into account
                        plot_npixels = 6.0
                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][0], pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][0]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][0])
                            x_range = (center_x - 0.5 * width, center_x + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0])
                            x_range = (center_x - width, center_x + width)
                        hist_residual_x_hist, hist_residual_x_xedges = np.histogram(difference[:, 0], range=x_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_x), np.max(intersection_x)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][1], npixels_per_bin * pixel_size[actual_dut][1])
                        else:
                            nbins = "auto"
                        _, hist_residual_x_yedges = np.histogram(intersection_x, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][1], pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][1]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][1]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][1])
                            y_range = (center_y - 0.5 * width, center_y + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_y / pixel_size[actual_dut][1])
                            y_range = (center_y - width, center_y + width)
                        hist_residual_y_hist, hist_residual_y_yedges = np.histogram(difference[:, 1], range=y_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_y), np.max(intersection_y)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][0], npixels_per_bin * pixel_size[actual_dut][0])
                        else:
                            nbins = "auto"
                        _, hist_residual_y_xedges = np.histogram(intersection_y, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][0], pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][0]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][0]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][0])
                            col_range = (center_col - 0.5 * width, center_col + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][0] * np.ceil(plot_npixels * fwhm_col / pixel_size[actual_dut][0])
                            col_range = (center_col - width, center_col + width)
                        hist_residual_col_hist, hist_residual_col_xedges = np.histogram(difference_local[:, 0], range=col_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_x_local), np.max(intersection_x_local)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][1], npixels_per_bin * pixel_size[actual_dut][1])
                        else:
                            nbins = "auto"
                        _, hist_residual_col_yedges = np.histogram(intersection_x_local, bins=nbins)

                        if nbins_per_pixel is not None:
                            width = max(plot_npixels * pixel_size[actual_dut][1], pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_x / pixel_size[actual_dut][1]))
                            if np.mod(width / pixel_size[actual_dut][0], 2) != 0:
                                width += pixel_size[actual_dut][1]
                            nbins = int(nbins_per_pixel * width / pixel_size[actual_dut][1])
                            row_range = (center_row - 0.5 * width, center_row + 0.5 * width)
                        else:
                            nbins = "auto"
                            width = pixel_size[actual_dut][1] * np.ceil(plot_npixels * fwhm_row / pixel_size[actual_dut][1])
                            row_range = (center_row - width, center_row + width)
                        hist_residual_row_hist, hist_residual_row_yedges = np.histogram(difference_local[:, 1], range=row_range, bins=nbins)

                        if npixels_per_bin is not None:
                            min_intersection, max_intersection = np.min(intersection_y_local), np.max(intersection_y_local)
                            nbins = np.arange(min_intersection, max_intersection + npixels_per_bin * pixel_size[actual_dut][0], npixels_per_bin * pixel_size[actual_dut][0])
                        else:
                            nbins = "auto"
                        _, hist_residual_row_xedges = np.histogram(intersection_y_local, bins=nbins)

                        # global x residual against x position
                        hist_x_residual_x_hist, hist_x_residual_x_xedges, hist_x_residual_x_yedges = np.histogram2d(
                            intersection_x,
                            difference[:, 0],
                            bins=(hist_residual_x_yedges, hist_residual_x_xedges))

                        # global y residual against y position
                        hist_y_residual_y_hist, hist_y_residual_y_xedges, hist_y_residual_y_yedges = np.histogram2d(
                            intersection_y,
                            difference[:, 1],
                            bins=(hist_residual_y_xedges, hist_residual_y_yedges))

                        # global y residual against x position
                        hist_x_residual_y_hist, hist_x_residual_y_xedges, hist_x_residual_y_yedges = np.histogram2d(
                            intersection_x,
                            difference[:, 1],
                            bins=(hist_residual_x_yedges, hist_residual_y_yedges))

                        # global x residual against y position
                        hist_y_residual_x_hist, hist_y_residual_x_xedges, hist_y_residual_x_yedges = np.histogram2d(
                            intersection_y,
                            difference[:, 0],
                            bins=(hist_residual_y_xedges, hist_residual_x_xedges))

                        # local column residual against column position
                        hist_col_residual_col_hist, hist_col_residual_col_xedges, hist_col_residual_col_yedges = np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 0],
                            bins=(hist_residual_col_yedges, hist_residual_col_xedges))

                        # local row residual against row position
                        hist_row_residual_row_hist, hist_row_residual_row_xedges, hist_row_residual_row_yedges = np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 1],
                            bins=(hist_residual_row_xedges, hist_residual_row_yedges))

                        # local row residual against column position
                        hist_col_residual_row_hist, hist_col_residual_row_xedges, hist_col_residual_row_yedges = np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 1],
                            bins=(hist_residual_col_yedges, hist_residual_row_yedges))

                        # local column residual against row position
                        hist_row_residual_col_hist, hist_row_residual_col_xedges, hist_row_residual_col_yedges = np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 0],
                            bins=(hist_residual_row_xedges, hist_residual_col_xedges))

                    else:  # adding data to existing histograms
                        hist_residual_x_hist += np.histogram(difference[:, 0], bins=hist_residual_x_xedges)[0]
                        hist_residual_y_hist += np.histogram(difference[:, 1], bins=hist_residual_y_yedges)[0]
                        hist_residual_col_hist += np.histogram(difference_local[:, 0], bins=hist_residual_col_xedges)[0]
                        hist_residual_row_hist += np.histogram(difference_local[:, 1], bins=hist_residual_row_yedges)[0]

                        # global x residual against x position
                        hist_x_residual_x_hist += np.histogram2d(
                            intersection_x,
                            difference[:, 0],
                            bins=(hist_x_residual_x_xedges, hist_x_residual_x_yedges))[0]

                        # global y residual against y position
                        hist_y_residual_y_hist += np.histogram2d(
                            intersection_y,
                            difference[:, 1],
                            bins=(hist_y_residual_y_xedges, hist_y_residual_y_yedges))[0]

                        # global y residual against x position
                        hist_x_residual_y_hist += np.histogram2d(
                            intersection_x,
                            difference[:, 1],
                            bins=(hist_x_residual_y_xedges, hist_x_residual_y_yedges))[0]

                        # global x residual against y position
                        hist_y_residual_x_hist += np.histogram2d(
                            intersection_y,
                            difference[:, 0],
                            bins=(hist_y_residual_x_xedges, hist_y_residual_x_yedges))[0]

                        # local column residual against column position
                        hist_col_residual_col_hist += np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 0],
                            bins=(hist_col_residual_col_xedges, hist_col_residual_col_yedges))[0]

                        # local row residual against row position
                        hist_row_residual_row_hist += np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 1],
                            bins=(hist_row_residual_row_xedges, hist_row_residual_row_yedges))[0]

                        # local row residual against column position
                        hist_col_residual_row_hist += np.histogram2d(
                            intersection_x_local,
                            difference_local[:, 1],
                            bins=(hist_col_residual_row_xedges, hist_col_residual_row_yedges))[0]

                        # local column residual against row position
                        hist_row_residual_col_hist += np.histogram2d(
                            intersection_y_local,
                            difference_local[:, 0],
                            bins=(hist_row_residual_col_xedges, hist_row_residual_col_yedges))[0]

                logging.debug('Storing residual histograms...')

                dut_name = dut_names[actual_dut] if dut_names else ("DUT" + str(actual_dut))
                # Global residuals
                fit_residual_x, cov_residual_x = analysis_utils.fit_residuals(
                    hist=hist_residual_x_hist,
                    edges=hist_residual_x_xedges,
                    label='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                      name='ResidualsX_DUT%d' % (actual_dut),
                                                      title='Residual distribution in x direction for %s' % (dut_name),
                                                      atom=tb.Atom.from_dtype(hist_residual_x_hist.dtype),
                                                      shape=hist_residual_x_hist.shape,
                                                      filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_x.attrs.xedges = hist_residual_x_xedges
                out_res_x.attrs.fit_coeff = fit_residual_x
                out_res_x.attrs.fit_cov = cov_residual_x
                out_res_x[:] = hist_residual_x_hist

                fit_residual_y, cov_residual_y = analysis_utils.fit_residuals(
                    hist=hist_residual_y_hist,
                    edges=hist_residual_y_yedges,
                    label='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                      name='ResidualsY_DUT%d' % (actual_dut),
                                                      title='Residual distribution in y direction for %s' % (dut_name),
                                                      atom=tb.Atom.from_dtype(hist_residual_y_hist.dtype),
                                                      shape=hist_residual_y_hist.shape,
                                                      filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_y.attrs.yedges = hist_residual_y_yedges
                out_res_y.attrs.fit_coeff = fit_residual_y
                out_res_y.attrs.fit_cov = cov_residual_y
                out_res_y[:] = hist_residual_y_hist

                fit_x_residual_x, cov_x_residual_x = analysis_utils.fit_residuals_vs_position(
                    hist=hist_x_residual_x_hist,
                    xedges=hist_x_residual_x_xedges,
                    yedges=hist_x_residual_x_yedges,
                    xlabel='X position [um]',
                    ylabel='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )

                out_x_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                        name='XResidualsX_DUT%d' % (actual_dut),
                                                        title='Residual distribution in x direction as a function of the x position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_x_residual_x_hist.dtype),
                                                        shape=hist_x_residual_x_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_x_res_x.attrs.xedges = hist_x_residual_x_xedges
                out_x_res_x.attrs.yedges = hist_x_residual_x_yedges
                out_x_res_x.attrs.fit_coeff = fit_x_residual_x
                out_x_res_x.attrs.fit_cov = cov_x_residual_x
                out_x_res_x[:] = hist_x_residual_x_hist

                fit_y_residual_y, cov_y_residual_y = analysis_utils.fit_residuals_vs_position(
                    hist=hist_y_residual_y_hist,
                    xedges=hist_y_residual_y_xedges,
                    yedges=hist_y_residual_y_yedges,
                    xlabel='Y position [um]',
                    ylabel='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_y_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                        name='YResidualsY_DUT%d' % (actual_dut),
                                                        title='Residual distribution in y direction as a function of the y position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_y_residual_y_hist.dtype),
                                                        shape=hist_y_residual_y_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_y_res_y.attrs.xedges = hist_y_residual_y_xedges
                out_y_res_y.attrs.yedges = hist_y_residual_y_yedges
                out_y_res_y.attrs.fit_coeff = fit_y_residual_y
                out_y_res_y.attrs.fit_cov = cov_y_residual_y
                out_y_res_y[:] = hist_y_residual_y_hist

                fit_x_residual_y, cov_x_residual_y = analysis_utils.fit_residuals_vs_position(
                    hist=hist_x_residual_y_hist,
                    xedges=hist_x_residual_y_xedges,
                    yedges=hist_x_residual_y_yedges,
                    xlabel='X position [um]',
                    ylabel='Y residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_x_res_y = out_file_h5.create_carray(out_file_h5.root,
                                                        name='XResidualsY_DUT%d' % (actual_dut),
                                                        title='Residual distribution in y direction as a function of the x position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_x_residual_y_hist.dtype),
                                                        shape=hist_x_residual_y_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_x_res_y.attrs.xedges = hist_x_residual_y_xedges
                out_x_res_y.attrs.yedges = hist_x_residual_y_yedges
                out_x_res_y.attrs.fit_coeff = fit_x_residual_y
                out_x_res_y.attrs.fit_cov = cov_x_residual_y
                out_x_res_y[:] = hist_x_residual_y_hist

                fit_y_residual_x, cov_y_residual_x = analysis_utils.fit_residuals_vs_position(
                    hist=hist_y_residual_x_hist,
                    xedges=hist_y_residual_x_xedges,
                    yedges=hist_y_residual_x_yedges,
                    xlabel='Y position [um]',
                    ylabel='X residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_y_res_x = out_file_h5.create_carray(out_file_h5.root,
                                                        name='YResidualsX_DUT%d' % (actual_dut),
                                                        title='Residual distribution in x direction as a function of the y position for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_y_residual_x_hist.dtype),
                                                        shape=hist_y_residual_x_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_y_res_x.attrs.xedges = hist_y_residual_x_xedges
                out_y_res_x.attrs.yedges = hist_y_residual_x_yedges
                out_y_res_x.attrs.fit_coeff = fit_y_residual_x
                out_y_res_x.attrs.fit_cov = cov_y_residual_x
                out_y_res_x[:] = hist_y_residual_x_hist

                # Local residuals
                fit_residual_col, cov_residual_col = analysis_utils.fit_residuals(
                    hist=hist_residual_col_hist,
                    edges=hist_residual_col_xedges,
                    label='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                        name='ResidualsCol_DUT%d' % (actual_dut),
                                                        title='Residual distribution in column direction for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_residual_col_hist.dtype),
                                                        shape=hist_residual_col_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_col.attrs.xedges = hist_residual_col_xedges
                out_res_col.attrs.fit_coeff = fit_residual_col
                out_res_col.attrs.fit_cov = cov_residual_col
                out_res_col[:] = hist_residual_col_hist

                fit_residual_row, cov_residual_row = analysis_utils.fit_residuals(
                    hist=hist_residual_row_hist,
                    edges=hist_residual_row_yedges,
                    label='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                        name='ResidualsRow_DUT%d' % (actual_dut),
                                                        title='Residual distribution in row direction for %s' % (dut_name),
                                                        atom=tb.Atom.from_dtype(hist_residual_row_hist.dtype),
                                                        shape=hist_residual_row_hist.shape,
                                                        filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_res_row.attrs.yedges = hist_residual_row_yedges
                out_res_row.attrs.fit_coeff = fit_residual_row
                out_res_row.attrs.fit_cov = cov_residual_row
                out_res_row[:] = hist_residual_row_hist

                fit_col_residual_col, cov_col_residual_col = analysis_utils.fit_residuals_vs_position(
                    hist=hist_col_residual_col_hist,
                    xedges=hist_col_residual_col_xedges,
                    yedges=hist_col_residual_col_yedges,
                    xlabel='Column position [um]',
                    ylabel='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_col_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                            name='ColResidualsCol_DUT%d' % (actual_dut),
                                                            title='Residual distribution in column direction as a function of the column position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_col_residual_col_hist.dtype),
                                                            shape=hist_col_residual_col_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_col_res_col.attrs.xedges = hist_col_residual_col_xedges
                out_col_res_col.attrs.yedges = hist_col_residual_col_yedges
                out_col_res_col.attrs.fit_coeff = fit_col_residual_col
                out_col_res_col.attrs.fit_cov = cov_col_residual_col
                out_col_res_col[:] = hist_col_residual_col_hist

                fit_row_residual_row, cov_row_residual_row = analysis_utils.fit_residuals_vs_position(
                    hist=hist_row_residual_row_hist,
                    xedges=hist_row_residual_row_xedges,
                    yedges=hist_row_residual_row_yedges,
                    xlabel='Row position [um]',
                    ylabel='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_row_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                            name='RowResidualsRow_DUT%d' % (actual_dut),
                                                            title='Residual distribution in row direction as a function of the row position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_row_residual_row_hist.dtype),
                                                            shape=hist_row_residual_row_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_row_res_row.attrs.xedges = hist_row_residual_row_xedges
                out_row_res_row.attrs.yedges = hist_row_residual_row_yedges
                out_row_res_row.attrs.fit_coeff = fit_row_residual_row
                out_row_res_row.attrs.fit_cov = cov_row_residual_row
                out_row_res_row[:] = hist_row_residual_row_hist

                fit_col_residual_row, cov_col_residual_row = analysis_utils.fit_residuals_vs_position(
                    hist=hist_col_residual_row_hist,
                    xedges=hist_col_residual_row_xedges,
                    yedges=hist_col_residual_row_yedges,
                    xlabel='Column position [um]',
                    ylabel='Row residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_col_res_row = out_file_h5.create_carray(out_file_h5.root,
                                                            name='ColResidualsRow_DUT%d' % (actual_dut),
                                                            title='Residual distribution in row direction as a function of the column position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_col_residual_row_hist.dtype),
                                                            shape=hist_col_residual_row_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_col_res_row.attrs.xedges = hist_col_residual_row_xedges
                out_col_res_row.attrs.yedges = hist_col_residual_row_yedges
                out_col_res_row.attrs.fit_coeff = fit_col_residual_row
                out_col_res_row.attrs.fit_cov = cov_col_residual_row
                out_col_res_row[:] = hist_col_residual_row_hist

                fit_row_residual_col, cov_row_residual_col = analysis_utils.fit_residuals_vs_position(
                    hist=hist_row_residual_col_hist,
                    xedges=hist_row_residual_col_xedges,
                    yedges=hist_row_residual_col_yedges,
                    xlabel='Row position [um]',
                    ylabel='Column residual [um]',
                    title='Residuals for %s' % (dut_name,),
                    output_pdf=output_pdf
                )
                out_row_res_col = out_file_h5.create_carray(out_file_h5.root,
                                                            name='RowResidualsCol_DUT%d' % (actual_dut),
                                                            title='Residual distribution in column direction as a function of the row position for %s' % (dut_name),
                                                            atom=tb.Atom.from_dtype(hist_row_residual_col_hist.dtype),
                                                            shape=hist_row_residual_col_hist.shape,
                                                            filters=tb.Filters(complib='blosc', complevel=5, fletcher32=False))
                out_row_res_col.attrs.xedges = hist_row_residual_col_xedges
                out_row_res_col.attrs.yedges = hist_row_residual_col_yedges
                out_row_res_col.attrs.fit_coeff = fit_row_residual_col
                out_row_res_col.attrs.fit_cov = cov_row_residual_col
                out_row_res_col[:] = hist_row_residual_col_hist

    if output_pdf is not None:
        output_pdf.close()


def calculate_efficiency(input_tracks_file, input_alignment_file, bin_size, sensor_size, output_efficiency_file=None, pixel_size=None, n_pixels=None, minimum_track_density=1, max_distance=500, use_duts=None, max_chi2=None, force_prealignment=False, cut_distance=None, col_range=None, row_range=None, show_inefficient_events=False, plot=True, chunk_size=1000000):
    '''Takes the tracks and calculates the hit efficiency and hit/track hit distance for selected DUTs.

    Parameters
    ----------
    input_tracks_file : string
        Filename of the input tracks file.
    input_alignment_file : string
        Filename of the input alignment file.
    bin_size : iterable
        Sizes of bins (i.e. (virtual) pixel size). Give one tuple (x, y) for every plane or list of tuples for different planes.
    sensor_size : Tuple or list of tuples
        Describes the sensor size for each DUT. If one tuple is given it is (size x, size y)
        If several tuples are given it is [(DUT0 size x, DUT0 size y), (DUT1 size x, DUT1 size y), ...]
    output_efficiency_file : string
        Filename of the output efficiency file. If None, the filename will be derived from the input hits file.
    minimum_track_density : int
        Minimum track density required to consider bin for efficiency calculation.
    use_duts : iterable
        Calculate the efficiency for selected DUTs. If None, all duts are selected.
    max_chi2 : uint
        Only use tracks with a chi2 <= max_chi2.
    force_prealignment : bool
        Take the prealignment, although if a coarse alignment is availale.
    cut_distance : int
        Use only distances (between DUT hit and track hit) smaller than cut_distance.
    max_distance : int
        Defines binnig of distance values.
    col_range, row_range : iterable
        Column / row value to calculate efficiency for (to neglect noisy edge pixels for efficiency calculation).
    plot : bool
        If True, create additional output plots.
    chunk_size : int
        Chunk size of the data when reading from file.
    '''
    logging.info('=== Calculating efficiency ===')

    if output_efficiency_file is None:
        output_efficiency_file = os.path.splitext(input_tracks_file)[0] + '_efficiency.h5'

    if plot is True:
        output_pdf = PdfPages(os.path.splitext(output_efficiency_file)[0] + '.pdf')
    else:
        output_pdf = None

    use_prealignment = True if force_prealignment else False

    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 not isinstance(max_chi2, Iterable):
        max_chi2 = [max_chi2] * len(use_duts)

    efficiencies = []
    pass_tracks = []
    total_tracks = []
    with tb.open_file(input_tracks_file, mode='r') as in_file_h5:
        for index, node in enumerate(in_file_h5.root):
            actual_dut = int(re.findall(r'\d+', node.name)[-1])
            if actual_dut not in use_duts:
                continue
            logging.info('Calculate efficiency for DUT%d', actual_dut)
{'fraction_non_alphanumeric': 0.06904313524475993, 'fraction_numerical': 0.007087061655207764, 'mean_word_length': 1.5093532421776696, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 348, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '805104', 'n_tokens_mistral': 25606, 'n_tokens_neox': 23923, 'n_words': 4587}