<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def fill_livetime_hist(skydir, tab_sc, tab_gti, zmax, costh_edges):
 """Generate a sequence of livetime distributions at the sky positions given by ``skydir``. The output of the method are two NxM arrays containing a sequence of histograms for N sky positions and M incidence angle bins where the bin edges are defined by ``costh_edges``. This method uses the same algorithm as `gtltcube` with the exception that SC time intervals are assumed to be aligned with GTIs. Parameters skydir : `~astropy.coordinates.SkyCoord` Vector of sky directions for which livetime histograms will be accumulated. tab_sc : `~astropy.table.Table` Spacecraft table. Must contain the following columns: START, STOP, LIVETIME, RA_SCZ, DEC_SZ, RA_ZENITH, DEC_ZENITH. tab_gti : `~astropy.table.Table` Table of good time intervals (GTIs). zmax : float Zenith cut. costh_edges : `~numpy.ndarray` Incidence angle bin edges in cos(angle). Returns ------- lt : `~numpy.ndarray` Array of livetime histograms. lt_wt : `~numpy.ndarray` Array of histograms of weighted livetime (livetime x livetime fraction). """


    if len(tab_gti) == 0:
        shape = (len(costh_edges) - 1, len(skydir))
        return (np.zeros(shape), np.zeros(shape))

    m = (tab_sc['START'] < tab_gti['STOP'][-1])
    m &= (tab_sc['STOP'] > tab_gti['START'][0])
    tab_sc = tab_sc[m]

    cos_zmax = np.cos(np.radians(zmax))
    sc_t0 = np.array(tab_sc['START'].data)
    sc_t1 = np.array(tab_sc['STOP'].data)
    sc_live = np.array(tab_sc['LIVETIME'].data)
    sc_lfrac = sc_live / (sc_t1 - sc_t0)

    sc_xyz = angle_to_cartesian(np.radians(tab_sc['RA_SCZ'].data),
                                np.radians(tab_sc['DEC_SCZ'].data))
    zn_xyz = angle_to_cartesian(np.radians(tab_sc['RA_ZENITH'].data),
                                np.radians(tab_sc['DEC_ZENITH'].data))

    tab_gti_t0 = np.array(tab_gti['START'].data)
    tab_gti_t1 = np.array(tab_gti['STOP'].data)

    # Index of the closest GTI interval
    idx = np.digitize(sc_t0, tab_gti_t0) - 1

    # start/stop time of closest GTI interval
    gti_t0 = np.zeros_like(sc_t0)
    gti_t1 = np.zeros_like(sc_t1)
    gti_t0[idx >= 0] = tab_gti_t0[idx[idx >= 0]]
    gti_t1[idx >= 0] = tab_gti_t1[idx[idx >= 0]]

    nbin = len(costh_edges) - 1

    lt = np.zeros((nbin,) + skydir.shape)
    lt_wt = np.zeros((nbin,) + skydir.shape)

    m0 = (idx >= 0) & (sc_t0 >= gti_t0) & (sc_t1 <= gti_t1)

    xyz = angle_to_cartesian(skydir.ra.rad, skydir.dec.rad)

    for i, t in enumerate(xyz):
        cos_sep = utils.dot_prod(t, sc_xyz)
        cos_zn = utils.dot_prod(t, zn_xyz)
        m = m0 & (cos_zn > cos_zmax) & (cos_sep > 0.0)
        bins = np.digitize(cos_sep[m], bins=costh_edges) - 1
        bins = np.clip(bins, 0, nbin - 1)
        lt[:, i] = np.bincount(bins, weights=sc_live[m], minlength=nbin)
        lt_wt[:, i] = np.bincount(bins, weights=sc_live[m] * sc_lfrac[m],
                                  minlength=nbin)

    return lt, lt_wt