Python astropy.units.TeV() Examples

def guess_shower_depth(energy):
    """
    Simple estimation of depth of shower max based on the expected gamma-ray elongation
    rate.

    Parameters
    ----------
    energy: float
        Energy of the shower in TeV

    Returns
    -------
    float: Expected depth of shower maximum
    """

    x_max_exp = 300 + 93 * np.log10(energy)

    return x_max_exp 

def predict_time(self, tel_type, energy, impact, x_max):
        """Creates predicted image for the specified pixels, interpolated
        from the template library.

        Parameters
        ----------
        tel_type: string
            Telescope type specifier
        energy: float
            Event energy (TeV)
        impact: float
            Impact diance of shower (metres)
        x_max: float
            Depth of shower maximum (num bins from expectation)

        Returns
        -------
        ndarray: predicted amplitude for all pixels

        """
        return self.time_prediction[tel_type](energy, impact, x_max) 

def test_write_container(temp_h5_file):
    r0tel = R0CameraContainer()
    mc = MCEventContainer()
    mc.reset()
    r0tel.waveform = np.random.uniform(size=(50, 10))
    r0tel.meta["test_attribute"] = 3.14159
    r0tel.meta["date"] = "2020-10-10"

    with HDF5TableWriter(
        temp_h5_file, group_name="R0", filters=tables.Filters(complevel=7)
    ) as writer:
        writer.exclude("tel_002", ".*samples")  # test exclusion of columns

        for ii in range(100):
            r0tel.waveform[:] = np.random.uniform(size=(50, 10))
            mc.energy = 10 ** np.random.uniform(1, 2) * u.TeV
            mc.core_x = np.random.uniform(-1, 1) * u.m
            mc.core_y = np.random.uniform(-1, 1) * u.m

            writer.write("tel_001", r0tel)
            writer.write("tel_002", r0tel)  # write a second table too
            writer.write("MC", mc) 

def __init__(
        self, regressor=RandomForestRegressor, cam_id_list="cam", unit=u.TeV, **kwargs
    ):
        super().__init__(model=regressor, cam_id_list=cam_id_list, unit=unit, **kwargs) 

def load(cls, path, cam_id_list, unit=u.TeV):
        """this is only here to overwrite the unit argument with an astropy
        quantity

        Parameters
        ----------
        path : string
            the path where the pre-trained, pickled regressors are
            stored `path` is assumed to contain a `{cam_id}` keyword
            to be replaced by each camera identifier in `cam_id_list`
            (or at least a naked `{}`).
        cam_id_list : list
            list of camera identifiers like telescope ID or camera ID
            and the assumed distinguishing feature in the filenames of
            the various pickled regressors.
        unit : astropy.Quantity
            scikit-learn regressor do not work with units. so append
            this one to the predictions. assuming that the models
            where trained with consistent units. (default: u.TeV)

        Returns
        -------
        EnergyRegressor:
            a ready-to-use instance of this class to predict any
            quantity you have trained for

        """
        return super().load(path, cam_id_list, unit) 

def image_prediction(self, tel_type, energy, impact, x_max, pix_x, pix_y):
        """Creates predicted image for the specified pixels, interpolated
        from the template library.

        Parameters
        ----------
        tel_type: string
            Telescope type specifier
        energy: float
            Event energy (TeV)
        impact: float
            Impact diance of shower (metres)
        x_max: float
            Depth of shower maximum (num bins from expectation)
        pix_x: ndarray
            X coordinate of pixels
        pix_y: ndarray
            Y coordinate of pixels

        Returns
        -------
        ndarray: predicted amplitude for all pixels

        """

        return self.prediction[tel_type](energy, impact, x_max, pix_x, pix_y) 

def test_showermaxestimator(en=5 * u.TeV, h_first_int=10 * u.km, az=70 * u.deg):
    estim = ShowerMaxEstimator(atmosphere_profile_name="paranal")
    h_max = estim.find_shower_max_height(en, h_first_int, az)
    assert h_max.unit.is_equivalent(u.m), "return value has not proper unit"
    return h_max 

def test_prepare_model():
    cam_id_list = ["FlashCam", "ASTRICam"]
    feature_list = {
        "FlashCam": [[1, 10], [2, 20], [3, 30], [0.9, 9],],
        "ASTRICam": [[10, 1], [20, 2], [30, 3], [9, 0.9],],
    }
    target_list = {
        "FlashCam": np.array([1, 2, 3, 0.9]) * u.TeV,
        "ASTRICam": np.array([1, 2, 3, 0.9]) * u.TeV,
    }

    reg = EnergyRegressor(cam_id_list=cam_id_list, n_estimators=10)
    reg.fit(feature_list, target_list)
    return reg, cam_id_list 

def _parse_mc_header(self):
        mc_run_head = self.file_.mc_run_headers[-1]

        return MCHeaderContainer(
            corsika_version=mc_run_head["shower_prog_vers"],
            simtel_version=mc_run_head["detector_prog_vers"],
            energy_range_min=mc_run_head["E_range"][0] * u.TeV,
            energy_range_max=mc_run_head["E_range"][1] * u.TeV,
            prod_site_B_total=mc_run_head["B_total"] * u.uT,
            prod_site_B_declination=Angle(mc_run_head["B_declination"], u.rad,),
            prod_site_B_inclination=Angle(mc_run_head["B_inclination"], u.rad,),
            prod_site_alt=mc_run_head["obsheight"] * u.m,
            spectral_index=mc_run_head["spectral_index"],
            shower_prog_start=mc_run_head["shower_prog_start"],
            shower_prog_id=mc_run_head["shower_prog_id"],
            detector_prog_start=mc_run_head["detector_prog_start"],
            detector_prog_id=mc_run_head["detector_prog_id"],
            num_showers=mc_run_head["n_showers"],
            shower_reuse=mc_run_head["n_use"],
            max_alt=mc_run_head["alt_range"][1] * u.rad,
            min_alt=mc_run_head["alt_range"][0] * u.rad,
            max_az=mc_run_head["az_range"][1] * u.rad,
            min_az=mc_run_head["az_range"][0] * u.rad,
            diffuse=mc_run_head["diffuse"],
            max_viewcone_radius=mc_run_head["viewcone"][1] * u.deg,
            min_viewcone_radius=mc_run_head["viewcone"][0] * u.deg,
            max_scatter_range=mc_run_head["core_range"][1] * u.m,
            min_scatter_range=mc_run_head["core_range"][0] * u.m,
            core_pos_mode=mc_run_head["core_pos_mode"],
            injection_height=mc_run_head["injection_height"] * u.m,
            atmosphere=mc_run_head["atmosphere"],
            corsika_iact_options=mc_run_head["corsika_iact_options"],
            corsika_low_E_model=mc_run_head["corsika_low_E_model"],
            corsika_high_E_model=mc_run_head["corsika_high_E_model"],
            corsika_bunchsize=mc_run_head["corsika_bunchsize"],
            corsika_wlen_min=mc_run_head["corsika_wlen_min"] * u.nm,
            corsika_wlen_max=mc_run_head["corsika_wlen_max"] * u.nm,
            corsika_low_E_detail=mc_run_head["corsika_low_E_detail"],
            corsika_high_E_detail=mc_run_head["corsika_high_E_detail"],
            run_array_direction=Angle(self.file_.header["direction"] * u.rad),
        ) 

def test_additional_meta_data_from_mc_header():
    with SimTelEventSource(input_url=gamma_test_large_path) as reader:
        data = next(iter(reader))

    # for expectation values
    from astropy import units as u
    from astropy.coordinates import Angle

    assert data.mcheader.corsika_version == 6990
    assert data.mcheader.spectral_index == -2.0
    assert data.mcheader.shower_reuse == 20
    assert data.mcheader.core_pos_mode == 1
    assert data.mcheader.diffuse == 1
    assert data.mcheader.atmosphere == 26

    # value read by hand from input card
    name_expectation = {
        "energy_range_min": u.Quantity(3.0e-03, u.TeV),
        "energy_range_max": u.Quantity(3.3e02, u.TeV),
        "prod_site_B_total": u.Quantity(23.11772346496582, u.uT),
        "prod_site_B_declination": Angle(0.0 * u.rad),
        "prod_site_B_inclination": Angle(-0.39641156792640686 * u.rad),
        "prod_site_alt": 2150.0 * u.m,
        "max_scatter_range": 3000.0 * u.m,
        "min_az": 0.0 * u.rad,
        "min_alt": 1.2217305 * u.rad,
        "max_viewcone_radius": 10.0 * u.deg,
        "corsika_wlen_min": 240 * u.nm,
    }

    for name, expectation in name_expectation.items():
        value = getattr(data.mcheader, name)

        assert value.unit == expectation.unit
        assert np.isclose(
            value.to_value(expectation.unit), expectation.to_value(expectation.unit)
        ) 

def evaluate(self, x, redshift):

            if isinstance(x, astropy_units.Quantity):

                # ebltable expects TeV
                eTeV = x.to(astropy_units.TeV).value
                return np.exp(-self._tau.opt_depth(redshift.value, eTeV)) * astropy_units.dimensionless_unscaled

            else:

                # otherwise it's in keV
                eTeV = old_div(x, 1e9)
                return np.exp(-self._tau.opt_depth(redshift, eTeV)) 

def test_fractional_powers():
    """See #2069"""
    m = 1e9 * u.Msun
    tH = 1. / (70. * u.km / u.s / u.Mpc)
    vc = 200 * u.km/u.s

    x = (c.G ** 2 * m ** 2 * tH.cgs) ** Fraction(1, 3) / vc
    v1 = x.to('pc')

    x = (c.G ** 2 * m ** 2 * tH) ** Fraction(1, 3) / vc
    v2 = x.to('pc')

    x = (c.G ** 2 * m ** 2 * tH.cgs) ** (1.0 / 3.0) / vc
    v3 = x.to('pc')

    x = (c.G ** 2 * m ** 2 * tH) ** (1.0 / 3.0) / vc
    v4 = x.to('pc')

    assert_allclose(v1, v2)
    assert_allclose(v2, v3)
    assert_allclose(v3, v4)

    x = u.m ** (1.0 / 101.0)
    assert isinstance(x.powers[0], float)

    x = u.m ** (3.0 / 7.0)
    assert isinstance(x.powers[0], Fraction)
    assert x.powers[0].numerator == 3
    assert x.powers[0].denominator == 7

    x = u.cm ** Fraction(1, 2) * u.cm ** Fraction(2, 3)
    assert isinstance(x.powers[0], Fraction)
    assert x.powers[0] == Fraction(7, 6)

    # Regression test for #9258.
    x = (u.TeV ** (-2.2)) ** (1/-2.2)
    assert isinstance(x.powers[0], Fraction)
    assert x.powers[0] == Fraction(1, 1) 

def test_image_prediction(self):
        pixel_x = np.array([0]) * u.deg
        pixel_y = np.array([0]) * u.deg

        image = np.array([1])
        pixel_area = np.array([1]) * u.deg * u.deg

        self.impact_reco.set_event_properties(
            {1: image},
            {1: pixel_x},
            {1: pixel_y},
            {1: pixel_area},
            {1: "CHEC"},
            {1: 0 * u.m},
            {1: 0 * u.m},
            array_direction=[0 * u.deg, 0 * u.deg],
        )

        """First check image prediction by directly accessing the function"""
        pred = self.impact_reco.image_prediction(
            "CHEC",
            zenith=0,
            azimuth=0,
            energy=1,
            impact=50,
            x_max=0,
            pix_x=pixel_x,
            pix_y=pixel_y,
        )

        assert np.sum(pred) != 0

        """Then check helper function gives the same answer"""
        shower = ReconstructedShowerContainer()
        shower.is_valid = True
        shower.alt = 0 * u.deg
        shower.az = 0 * u.deg
        shower.core_x = 0 * u.m
        shower.core_y = 100 * u.m
        shower.h_max = 300 + 93 * np.log10(1)

        energy = ReconstructedEnergyContainer()
        energy.is_valid = True
        energy.energy = 1 * u.TeV
        pred2 = self.impact_reco.get_prediction(
            1, shower_reco=shower, energy_reco=energy
        )
        print(pred, pred2)
        assert pred.all() == pred2.all() 

def _write_simulation_histograms(self, writer: HDF5TableWriter):
        """ Write the distribution of thrown showers

        Notes
        -----
        - this only runs if this is a simulation file. The current implementation is
          a bit of a hack and implies we should improve SimTelEventSource to read this
          info.
        - Currently the histograms are at the end of the simtel file, so if max_events
          is set to non-zero, the end of the file may not be read, and this no
          histograms will be found.
        """
        self.log.debug("Writing simulation histograms")

        if not isinstance(self.event_source, SimTelEventSource):
            return

        def fill_from_simtel(
            obs_id, eventio_hist, container: SimulatedShowerDistribution
        ):
            """ fill from a SimTel Histogram entry"""
            container.obs_id = obs_id
            container.hist_id = eventio_hist["id"]
            container.num_entries = eventio_hist["entries"]
            xbins = np.linspace(
                eventio_hist["lower_x"],
                eventio_hist["upper_x"],
                eventio_hist["n_bins_x"] + 1,
            )
            ybins = np.linspace(
                eventio_hist["lower_y"],
                eventio_hist["upper_y"],
                eventio_hist["n_bins_y"] + 1,
            )

            container.bins_core_dist = xbins * u.m
            container.bins_energy = 10 ** ybins * u.TeV
            container.histogram = eventio_hist["data"]
            container.meta["hist_title"] = eventio_hist["title"]
            container.meta["x_label"] = "Log10 E (TeV)"
            container.meta["y_label"] = "3D Core Distance (m)"

        hists = self.event_source.file_.histograms
        if hists is not None:
            hist_container = SimulatedShowerDistribution()
            hist_container.prefix = ""
            for hist in hists:
                if hist["id"] == 6:
                    fill_from_simtel(self.event_source.obs_id, hist, hist_container)
                    writer.write(
                        table_name="simulation/service/shower_distribution",
                        containers=hist_container,
                    ) 

def test_one_free_parameter_input_output():

    fluxUnit = 1.0 / (u.TeV * u.cm ** 2 * u.s)

    temp_file = "__test_mle.fits"

    spectrum = Powerlaw()
    source = PointSource("tst", ra=100, dec=20, spectral_shape=spectrum)
    model = Model(source)

    spectrum.piv = 7 * u.TeV
    spectrum.index = -2.3
    spectrum.K = 1e-15 * fluxUnit

    spectrum.piv.fix = True

    # two free parameters (one with units)
    spectrum.index.fix = False
    spectrum.K.fix = False
    cov_matrix = np.diag([0.001] * 2)
    ar = MLEResults(model, cov_matrix, {})

    ar.write_to(temp_file, overwrite=True)
    ar_reloaded = load_analysis_results(temp_file)
    os.remove(temp_file)
    _results_are_same(ar, ar_reloaded)

    # one free parameter with units
    spectrum.index.fix = True
    spectrum.K.fix = False
    cov_matrix = np.diag([0.001] * 1)
    ar = MLEResults(model, cov_matrix, {})

    ar.write_to(temp_file, overwrite=True)
    ar_reloaded = load_analysis_results(temp_file)
    os.remove(temp_file)
    _results_are_same(ar, ar_reloaded)

    # one free parameter without units
    spectrum.index.fix = False
    spectrum.K.fix = True
    cov_matrix = np.diag([0.001] * 1)
    ar = MLEResults(model, cov_matrix, {})

    ar.write_to(temp_file, overwrite=True)
    ar_reloaded = load_analysis_results(temp_file)
    os.remove(temp_file)
    _results_are_same(ar, ar_reloaded)