How to use the astropy.units.Quantity function in astropy

def test_fill_cube():
    filename = gammapy_extra.filename('test_datasets/background/bg_cube_model_test1.fits')
    array = Cube.read(filename, format='table', scheme='bg_cube')
    array.data = Quantity(np.zeros_like(array.data.value), 'u')
    print(type(array.data))

    dir = str(gammapy_extra.dir) + '/datasets/hess-crab4-hd-hap-prod2'
    data_store = DataStore.from_dir(dir)
    ev_list = data_store.load_all('events')

    array.fill_events(ev_list)

    array.write('test_background.fits', format='image', clobber=True)

The beginning of the Sedov-Taylor phase of the SNR is defined by the condition,
        that the swept up mass of the surrounding medium equals the mass of the
        ejected mass.

        The time scale is given by:

        .. math::
            t_{begin} \approx 200
            \left(\frac{E_{SN}}{10^{51}erg}\right)^{-1/2}
            \left(\frac{M_{ej}}{M_{\odot}}\right)^{5/6}
            \left(\frac{\rho_{ISM}}{10^{-24}g/cm^3}\right)^{-1/3}
            \text{yr}
        """
        term1 = (self.e_sn / Quantity(1e51, "erg")) ** (-1.0 / 2)
        term2 = (self.m_ejecta / const.M_sun) ** (5.0 / 6)
        term3 = (self.rho_ISM / (Quantity(1, "cm-3") * const.m_p)) ** (-1.0 / 3)
        return Quantity(200, "yr") * term1 * term2 * term3

# determine clock rate
            if adc_name == 'SIS 3305':
                # has different clock rate modes
                try:
                    cr_mode = config_group[name].attrs['Channel mode']
                    cr_mode = int(cr_mode)
                except (KeyError, ValueError):
                    why = ("HDF5 structure unexpected..."
                           + "'{}'".format(config_name + '/' + name)
                           + " does not define a clock rate mode"
                           + "...setting to None in the `configs` dict")
                    warn(why)
                    cr_mode = -1
                if cr_mode == 0:
                    cr = u.Quantity(1.25, unit='GHz')
                elif cr_mode == 1:
                    cr = u.Quantity(2.5, unit='GHz')
                elif cr_mode == 2:
                    cr = u.Quantity(5.0, unit='GHz')
                else:
                    cr = None
            else:
                # 'SIS 3302' has one clock rate mode
                cr = u.Quantity(100.0, unit='MHz')

            # build subconn tuple with connected board, channels, and
            # acquisition parameters
            subconn = (brd, tuple(chs),
                       {'clock rate': cr,
                        'shot average (software)': shot_ave,
                        'sample average (hardware)': sample_ave})

# Load exoplanets table
        table = ExoplanetOrbitDatabase.get_table(cache=cache, show_progress=show_progress)

        if not exoplanet_name.lower().strip() in table['NAME_LOWERCASE'].data:
            raise ValueError('Planet "{0}" not found in exoplanets.org catalog')

        row = table.loc[exoplanet_name.lower().strip()]

        kwargs = dict(parameter_source="exoplanets.org")
        for column in row.colnames:
            value = row[column]

            # If param is unitful quantity, make it a `astropy.units.Quantity`
            if table[column].unit is not None:
                parameter = u.Quantity(value, unit=table[column].unit)

            # If param describes a time, make it a `astropy.time.Time`
            elif column in TIME_ATTRS:
                parameter = Time(value, format=TIME_ATTRS[column])

            elif column in BOOL_ATTRS:
                parameter = bool(value)

            # Otherwise, simply set the parameter to its raw value
            else:
                parameter = value

            # Attributes should be all lowercase:
            kwargs[column.lower()] = parameter

        return cls(**kwargs)

# Defaults and input handling
        if theta is None:
            theta = Angle(0, "deg")
        else:
            theta = Angle(theta)

        if rad is None:
            rad = Angle(np.arange(0, 1.5, 0.005), "deg")
        else:
            rad = Angle(rad).to("deg")

        psf_value = Quantity(np.zeros((energies.size, rad.size)), "deg^-2")

        for idx, energy in enumerate(energies):
            psf_gauss = self.psf_at_energy_and_theta(energy, theta)
            psf_value[idx] = Quantity(psf_gauss(rad), "deg^-2")

        return EnergyDependentTablePSF(
            energy=energies, rad=rad, exposure=exposure, psf_value=psf_value
        )

pix_type, pix_rot = CameraGeometry.simtel_shape_to_type(pixel_shape)
        except ValueError:
            warnings.warn(
                f'Unkown pixel_shape {pixel_shape} for tel_id {tel_id}',
                UnknownPixelShapeWarning,
            )
            pix_type = 'hexagon'
            pix_rot = '0d'

        pix_area = u.Quantity(file.get_pixel_area(tel_id), u.m**2)

        mirror_area = u.Quantity(file.get_mirror_area(tel_id), u.m**2)
        num_tiles = file.get_mirror_number(tel_id)
        cam_rot = file.get_camera_rotation_angle(tel_id)
        num_mirrors = file.get_mirror_number(tel_id)
        sampling_rate = u.Quantity(1 / file.get_time_slice(tel_id), u.GHz)
        reference_pulse_shape = file.get_ref_shapes(tel_id)
        reference_pulse_sample_width = u.Quantity(file.get_ref_step(tel_id), u.ns)

        geometry = CameraGeometry(
            telescope.camera_name,
            pix_id=np.arange(n_pixels),
            pix_x=pix_x,
            pix_y=pix_y,
            pix_area=pix_area,
            pix_type=pix_type,
            pix_rotation=pix_rot,
            cam_rotation=-Angle(cam_rot, u.rad),
            apply_derotation=True,
        )
        readout = CameraReadout(
            telescope.camera_name,

cutout = np.copy(data[self.bbox.slices])
            else:
                cutout = data[self.bbox.slices]

        if partial_overlap or (cutout.shape != self.shape):
            slices_large, slices_small = self._overlap_slices(data.shape)

            if slices_small is None:
                return None    # no overlap

            # cutout is a copy
            cutout = np.zeros(self.shape, dtype=data.dtype)
            cutout[:] = fill_value
            cutout[slices_small] = data[slices_large]

            if isinstance(data, u.Quantity):
                cutout = u.Quantity(cutout, unit=data.unit)

        return cutout

-----------
        hillas_dict: dict[HillasContainer]
            dictionary of hillas moments
        array_pointing: SkyCoord[HorizonFrame]
            Pointing direction of the array

        Returns
        -----------
        core_x: u.Quantity
            estimated x position of impact
        core_y: u.Quantity
            estimated y position of impact

        """
        if self.divergent_mode:
            psi = u.Quantity(list(self.corrected_angle_dict.values()))
        else:
            psi = u.Quantity([h.psi for h in hillas_dict.values()])

        z = np.zeros(len(psi))
        uvw_vectors = np.column_stack([np.cos(psi).value, np.sin(psi).value, z])

        tilted_frame = TiltedGroundFrame(pointing_direction=array_pointing)
        ground_frame = GroundFrame()

        positions = [
            (
                SkyCoord(*plane.pos, frame=ground_frame)
                .transform_to(tilted_frame)
                .cartesian.xyz
            )
            for plane in self.hillas_planes.values()

rad : `~astropy.coordinates.Angle`
            Offset wrt source position

        Returns
        -------
        values : `~astropy.units.Quantity`
            Interpolated value
        """
        if energy is None:
            energy = self._energy_logcenter()
        if offset is None:
            offset = self.offset
        if rad is None:
            rad = self._rad_center()

        rad = np.atleast_1d(u.Quantity(rad))
        offset = np.atleast_1d(u.Quantity(offset))
        energy = np.atleast_1d(u.Quantity(energy))
        return self._interpolate(
            (
                rad[:, np.newaxis, np.newaxis],
                offset[np.newaxis, :, np.newaxis],
                energy[np.newaxis, np.newaxis, :],
            )

else:
            freq = self.periodogram.frequency  # Makes code below more readable
            power = self.periodogram.power     # Makes code below more readable

        fmin = freq[0]
        fmax = freq[-1]

        # Check for any superfluous input
        if (numax is not None) & (any([a is not None for a in [minimum_frequency, maximum_frequency]])):
            warnings.warn("You have passed both a numax and a frequency limit. "
                          "The frequency limit will override the numax input.",
                          LightkurveWarning)

        # Ensure input numax is in the correct units (if there is one)
        if numax is not None:
            numax = u.Quantity(numax, freq.unit).value
            if numax > freq[-1].value:
                raise ValueError("You can't pass in a numax outside the"
                                "frequency range of the periodogram.")

            fwhm = utils.get_fwhm(self.periodogram, numax)

            fmin = numax - 2*fwhm
            if fmin < freq[0].value:
                fmin = freq[0].value

            fmax = numax + 2*fwhm
            if fmax > freq[-1].value:
                fmax = freq[-1].value

        # Set limits and set them in the right units
        if minimum_frequency is not None: