pixel_clusterizer/cluster_functions.py
def _cluster_hits(hits, clusters, assigned_hit_array, cluster_hit_indices, column_cluster_distance, row_cluster_distance, frame_cluster_distance, min_hit_charge, max_hit_charge, ignore_same_hits, noisy_pixels, disabled_pixels):
    ''' Main precompiled function that loopes over the hits and clusters them
    '''
    total_hits = hits.shape[0]
    if total_hits == 0:
        return 0  # total clusters
    max_cluster_hits = cluster_hit_indices.shape[0]

    if total_hits != clusters.shape[0]:
        raise ValueError("hits and clusters must be the same size")

    if total_hits != assigned_hit_array.shape[0]:
        raise ValueError("hits and assigned_hit_array must be the same size")

    # Correction for charge weighting
    # Some chips have non-zero charge for a charge value of zero, charge needs to be corrected to calculate cluster center correctly
    if min_hit_charge == 0:
        charge_correction = 1
    else:
        charge_correction = 0

    # Temporary variables that are reset for each cluster or event
    start_event_hit_index = 0
    start_event_cluster_index = 0
    cluster_size = 0
    event_number = hits[0]['event_number']
    event_cluster_index = 0

    # Outer loop over all hits in the array (referred to as actual hit)
    for i in range(total_hits):

        # Check for new event and reset event variables
        if _new_event(hits[i]['event_number'], event_number):
            _finish_event(
                hits=hits,
                clusters=clusters,
                start_event_hit_index=start_event_hit_index,
                stop_event_hit_index=i,
                start_event_cluster_index=start_event_cluster_index,
                stop_event_cluster_index=start_event_cluster_index + event_cluster_index)

            start_event_hit_index = i
            start_event_cluster_index = start_event_cluster_index + event_cluster_index
            event_number = hits[i]['event_number']
            event_cluster_index = 0

        if assigned_hit_array[i] > 0:  # Hit was already assigned to a cluster in the inner loop, thus skip actual hit
            continue

        if not _hit_ok(
                hit=hits[i],
                min_hit_charge=min_hit_charge,
                max_hit_charge=max_hit_charge) or (disabled_pixels.shape[0] != 0 and _pixel_masked(hits[i], disabled_pixels)):
            _set_hit_invalid(hit=hits[i], cluster_id=-1)
            assigned_hit_array[i] = 1
            continue

        # Set/reset cluster variables for new cluster
        # Reset temp array with hit indices of actual cluster for the next cluster
        _set_1d_array(cluster_hit_indices, -1, cluster_size)
        cluster_hit_indices[0] = i
        assigned_hit_array[i] = 1
        cluster_size = 1  # actual cluster has one hit so far

        for j in cluster_hit_indices:  # Loop over all hits of the actual cluster; cluster_hit_indices is updated within the loop if new hit are found
            if j < 0:  # There are no more cluster hits found
                break

            for k in range(cluster_hit_indices[0] + 1, total_hits):
                # Stop event hits loop if new event is reached
                if _new_event(hits[k]['event_number'], event_number):
                    break

                # Hit is already assigned to a cluster, thus skip actual hit
                if assigned_hit_array[k] > 0:
                    continue

                if not _hit_ok(
                        hit=hits[k],
                        min_hit_charge=min_hit_charge,
                        max_hit_charge=max_hit_charge) or (disabled_pixels.shape[0] != 0 and _pixel_masked(hits[k], disabled_pixels)):
                    _set_hit_invalid(hit=hits[k], cluster_id=-1)
                    assigned_hit_array[k] = 1
                    continue

                # Check if event hit belongs to actual hit and thus to the actual cluster
                if _is_in_max_difference(hits[j]['column'], hits[k]['column'], column_cluster_distance) and _is_in_max_difference(hits[j]['row'], hits[k]['row'], row_cluster_distance) and _is_in_max_difference(hits[j]['frame'], hits[k]['frame'], frame_cluster_distance):
                    if not ignore_same_hits or hits[j]['column'] != hits[k]['column'] or hits[j]['row'] != hits[k]['row']:
                        cluster_size += 1
                        if cluster_size > max_cluster_hits:
                            raise IndexError('cluster_hit_indices is too small to contain all cluster hits')
                        cluster_hit_indices[cluster_size - 1] = k
                        assigned_hit_array[k] = 1

                    else:
                        _set_hit_invalid(hit=hits[k], cluster_id=-2)
                        assigned_hit_array[k] = 1

        # check for valid cluster and add it to the array
        if cluster_size == 1 and noisy_pixels.shape[0] != 0 and _pixel_masked(hits[cluster_hit_indices[0]], noisy_pixels):
            _set_hit_invalid(hit=hits[cluster_hit_indices[0]], cluster_id=-1)
        else:
            _finish_cluster(
                hits=hits,
                clusters=clusters,
                cluster_size=cluster_size,
                cluster_hit_indices=cluster_hit_indices,
                cluster_index=start_event_cluster_index + event_cluster_index,
                cluster_id=event_cluster_index,
                charge_correction=charge_correction,
                noisy_pixels=noisy_pixels,
                disabled_pixels=disabled_pixels)
            event_cluster_index += 1

    # Last event is assumed to be finished at the end of the hit array, thus add info
    _finish_event(
        hits=hits,
        clusters=clusters,
        start_event_hit_index=start_event_hit_index,
        stop_event_hit_index=total_hits,
        start_event_cluster_index=start_event_cluster_index,
        stop_event_cluster_index=start_event_cluster_index + event_cluster_index)

    total_clusters = start_event_cluster_index + event_cluster_index
    return total_clusters
def _cluster_hits(hits, clusters, assigned_hit_array, cluster_hit_indices, column_cluster_distance, row_cluster_distance, frame_cluster_distance, min_hit_charge, max_hit_charge, ignore_same_hits, noisy_pixels, disabled_pixels):
    ''' Main precompiled function that loopes over the hits and clusters them
    '''
    total_hits = hits.shape[0]
    if total_hits == 0:
        return 0  # total clusters
    max_cluster_hits = cluster_hit_indices.shape[0]

    if total_hits != clusters.shape[0]:
        raise ValueError("hits and clusters must be the same size")

    if total_hits != assigned_hit_array.shape[0]:
        raise ValueError("hits and assigned_hit_array must be the same size")

    # Correction for charge weighting
    # Some chips have non-zero charge for a charge value of zero, charge needs to be corrected to calculate cluster center correctly
    if min_hit_charge == 0:
        charge_correction = 1
    else:
        charge_correction = 0

    # Temporary variables that are reset for each cluster or event
    start_event_hit_index = 0
    start_event_cluster_index = 0
    cluster_size = 0
    event_number = hits[0]['event_number']
    event_cluster_index = 0

    # Outer loop over all hits in the array (referred to as actual hit)
    for i in range(total_hits):

        # Check for new event and reset event variables
        if _new_event(hits[i]['event_number'], event_number):
            _finish_event(
                hits=hits,
                clusters=clusters,
                start_event_hit_index=start_event_hit_index,
                stop_event_hit_index=i,
                start_event_cluster_index=start_event_cluster_index,
                stop_event_cluster_index=start_event_cluster_index + event_cluster_index)

            start_event_hit_index = i
            start_event_cluster_index = start_event_cluster_index + event_cluster_index
            event_number = hits[i]['event_number']
            event_cluster_index = 0

        if assigned_hit_array[i] > 0:  # Hit was already assigned to a cluster in the inner loop, thus skip actual hit
            continue

        if not _hit_ok(
                hit=hits[i],
                min_hit_charge=min_hit_charge,
                max_hit_charge=max_hit_charge) or (disabled_pixels.shape[0] != 0 and _pixel_masked(hits[i], disabled_pixels)):
            _set_hit_invalid(hit=hits[i], cluster_id=-1)
            assigned_hit_array[i] = 1
            continue

        # Set/reset cluster variables for new cluster
        # Reset temp array with hit indices of actual cluster for the next cluster
        _set_1d_array(cluster_hit_indices, -1, cluster_size)
        cluster_hit_indices[0] = i
        assigned_hit_array[i] = 1
        cluster_size = 1  # actual cluster has one hit so far

        for j in cluster_hit_indices:  # Loop over all hits of the actual cluster; cluster_hit_indices is updated within the loop if new hit are found
            if j < 0:  # There are no more cluster hits found
                break

            for k in range(cluster_hit_indices[0] + 1, total_hits):
                # Stop event hits loop if new event is reached
                if _new_event(hits[k]['event_number'], event_number):
                    break

                # Hit is already assigned to a cluster, thus skip actual hit
                if assigned_hit_array[k] > 0:
                    continue

                if not _hit_ok(
                        hit=hits[k],
                        min_hit_charge=min_hit_charge,
                        max_hit_charge=max_hit_charge) or (disabled_pixels.shape[0] != 0 and _pixel_masked(hits[k], disabled_pixels)):
                    _set_hit_invalid(hit=hits[k], cluster_id=-1)
                    assigned_hit_array[k] = 1
                    continue

                # Check if event hit belongs to actual hit and thus to the actual cluster
                if _is_in_max_difference(hits[j]['column'], hits[k]['column'], column_cluster_distance) and _is_in_max_difference(hits[j]['row'], hits[k]['row'], row_cluster_distance) and _is_in_max_difference(hits[j]['frame'], hits[k]['frame'], frame_cluster_distance):
                    if not ignore_same_hits or hits[j]['column'] != hits[k]['column'] or hits[j]['row'] != hits[k]['row']:
                        cluster_size += 1
                        if cluster_size > max_cluster_hits:
                            raise IndexError('cluster_hit_indices is too small to contain all cluster hits')
                        cluster_hit_indices[cluster_size - 1] = k
                        assigned_hit_array[k] = 1

                    else:
                        _set_hit_invalid(hit=hits[k], cluster_id=-2)
                        assigned_hit_array[k] = 1

        # check for valid cluster and add it to the array
        if cluster_size == 1 and noisy_pixels.shape[0] != 0 and _pixel_masked(hits[cluster_hit_indices[0]], noisy_pixels):
            _set_hit_invalid(hit=hits[cluster_hit_indices[0]], cluster_id=-1)
        else:
            _finish_cluster(
                hits=hits,
                clusters=clusters,
                cluster_size=cluster_size,
                cluster_hit_indices=cluster_hit_indices,
                cluster_index=start_event_cluster_index + event_cluster_index,
                cluster_id=event_cluster_index,
                charge_correction=charge_correction,
                noisy_pixels=noisy_pixels,
                disabled_pixels=disabled_pixels)
            event_cluster_index += 1

    # Last event is assumed to be finished at the end of the hit array, thus add info
    _finish_event(
        hits=hits,
        clusters=clusters,
        start_event_hit_index=start_event_hit_index,
        stop_event_hit_index=total_hits,
        start_event_cluster_index=start_event_cluster_index,
        stop_event_cluster_index=start_event_cluster_index + event_cluster_index)

    total_clusters = start_event_cluster_index + event_cluster_index
    return total_clusters
Main precompiled function that loopes over the hits and clusters them