Python GraphArtist.set_xlimits方法代码示例

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_trace(station_group, idx):
    events = station_group.events
    blobs = station_group.blobs

    traces_idx = events[idx]['traces']
    traces = get_traces(blobs, traces_idx)
    traces = array(traces)
    x = arange(traces.shape[1])
    x *= 2.5

    clf()
    plot(x, traces.T)
    xlim(0, 200)

    #line_styles = ['solid', 'dashed', 'dotted', 'dashdotted']
    line_styles = ['black', 'black!80', 'black!60', 'black!40']
    styles = (u for u in line_styles)

    graph = GraphArtist(width=r'.5\linewidth')
    for trace in traces:
        graph.plot(x, trace / 1000, mark=None, linestyle=styles.next())
    graph.set_xlabel(r"Time [\si{\nano\second}]")
    graph.set_ylabel(r"Signal [\si{\volt}]")
    graph.set_xlimits(0, 200)
    graph.save('plots/traces')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_arrival_times():
    graph = GraphArtist()

    figure()
    sim = data.root.showers.E_1PeV.zenith_22_5
    t = get_front_arrival_time(sim, 20, 5, pi / 8)
    n, bins = histogram(t, bins=linspace(0, 50, 201))
    mct = monte_carlo_timings(n, bins, 100000)
    n, bins, patches = hist(mct, bins=linspace(0, 20, 101), histtype='step')
    graph.histogram(n, bins, linestyle='black!50')

    mint = my_t_draw_something(data, 2, 100000)
    n, bins, patches = hist(mint, bins=linspace(0, 20, 101), histtype='step')
    graph.histogram(n, bins)

    xlabel("Arrival time [ns]")
    ylabel("Number of events")

    graph.set_xlabel(r"Arrival time [\si{\nano\second}]")
    graph.set_ylabel("Number of events")
    graph.set_xlimits(0, 20)
    graph.set_ylimits(min=0)
    graph.save('plots/SIM-T')

    print(median(t), median(mct), median(mint))

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def boxplot_arrival_times(group, N):
    table = group.E_1PeV.zenith_0

    sel = table.read_where('min_n134 >= N')
    t1 = sel[:]['t1']
    t3 = sel[:]['t3']
    t4 = sel[:]['t4']
    ts = concatenate([t1, t3, t4])
    print "Median arrival time delay over all detected events", median(ts)

    figure()

    bin_edges = linspace(0, 100, 11)
    x, arrival_times = [], []
    t25, t50, t75 = [], [], []
    for low, high in zip(bin_edges[:-1], bin_edges[1:]):
        query = '(min_n134 >= N) & (low <= r) & (r < high)'
        sel = table.read_where(query)
        t1 = sel[:]['t1']
        t2 = sel[:]['t2']
        ct1 = t1.compress((t1 > -999) & (t2 > -999))
        ct2 = t2.compress((t1 > -999) & (t2 > -999))
        ts = abs(ct2 - ct1)

        t25.append(scoreatpercentile(ts, 25))
        t50.append(scoreatpercentile(ts, 50))
        t75.append(scoreatpercentile(ts, 75))
        x.append((low + high) / 2)

    fill_between(x, t25, t75, color='0.75')
    plot(x, t50, 'o-', color='black')

    xlabel("Core distance [m]")
    ylabel("Arrival time delay [ns]")
    #title(r"$N_{MIP} \geq %d, \quad \theta = 0^\circ$" % N)

    xticks(arange(0, 100.5, 10))

    utils.savedata((x, t25, t50, t75), N)
    utils.saveplot(N)

    graph = GraphArtist()
    graph.shade_region(x, t25, t75)
    graph.plot(x, t50, linestyle=None)
    graph.set_xlabel(r"Core distance [\si{\meter}]")
    graph.set_ylabel(r"Arrival time difference $|t_2 - t_1|$ [\si{\nano\second}]")
    graph.set_xlimits(0, 100)
    graph.set_ylimits(min=0)
    artist.utils.save_graph(graph, suffix=N, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def boxplot_phi_reconstruction_results_for_MIP(table, N):
    figure()

    THETA = deg2rad(22.5)
    DTHETA = deg2rad(5.)

    bin_edges = linspace(-180, 180, 18)
    x, r_dphi = [], []
    d25, d50, d75 = [], [], []
    for low, high in zip(bin_edges[:-1], bin_edges[1:]):
        rad_low = deg2rad(low)
        rad_high = deg2rad(high)
        query = '(min_n134 >= N) & (rad_low < reference_phi) & (reference_phi < rad_high) & (abs(reference_theta - THETA) <= DTHETA)'
        sel = table.read_where(query)
        dphi = sel[:]['reconstructed_phi'] - sel[:]['reference_phi']
        dphi = (dphi + pi) % (2 * pi) - pi
        r_dphi.append(rad2deg(dphi))

        d25.append(scoreatpercentile(rad2deg(dphi), 25))
        d50.append(scoreatpercentile(rad2deg(dphi), 50))
        d75.append(scoreatpercentile(rad2deg(dphi), 75))
        x.append((low + high) / 2)

    #boxplot(r_dphi, positions=x, widths=1 * (high - low), sym='')
    fill_between(x, d25, d75, color='0.75')
    plot(x, d50, 'o-', color='black')

    xlabel(r"$\phi_K$ [deg]")
    ylabel(r"$\phi_H - \phi_K$ [deg]")
    title(r"$N_{MIP} \geq %d, \quad \theta = 22.5^\circ \pm %d^\circ$" % (N, rad2deg(DTHETA)))

    xticks(linspace(-180, 180, 9))
    axhline(0, color='black')

    utils.saveplot(N)

    graph = GraphArtist()
    graph.draw_horizontal_line(0, linestyle='gray')
    graph.shade_region(x, d25, d75)
    graph.plot(x, d50, linestyle=None)
    graph.set_xlabel(r"$\phi_K$ [\si{\degree}]")
    graph.set_ylabel(r"$\phi_H - \phi_K$ [\si{\degree}]")
    graph.set_xticks([-180, -90, '...', 180])
    graph.set_xlimits(-180, 180)
    graph.set_ylimits(-23, 23)
    artist.utils.save_graph(graph, suffix=N, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_detection_efficiency_vs_R_for_angles(N):
    figure()
    graph = GraphArtist()
    locations = iter(['right', 'left', 'below left'])
    positions = iter([.18, .14, .15])

    bin_edges = linspace(0, 100, 20)
    x = (bin_edges[:-1] + bin_edges[1:]) / 2.

    for angle in [0, 22.5, 35]:
        angle_str = str(angle).replace('.', '_')
        shower_group = '/simulations/E_1PeV/zenith_%s' % angle_str

        efficiencies = []
        for low, high in zip(bin_edges[:-1], bin_edges[1:]):
            shower_results = []
            for shower in data.list_nodes(shower_group):
                sel_query = '(low <= r) & (r < high)'
                coinc_sel = shower.coincidences.read_where(sel_query)
                ids = coinc_sel['id']
                obs_sel = shower.observables.read_coordinates(ids)
                assert (obs_sel['id'] == ids).all()

                o = obs_sel
                sel = obs_sel.compress((o['n1'] >= N) & (o['n3'] >= N) &
                                       (o['n4'] >= N))
                shower_results.append(len(sel) / len(obs_sel))
            efficiencies.append(mean(shower_results))

        plot(x, efficiencies, label=r'$\theta = %s^\circ$' % angle)
        graph.plot(x, efficiencies, mark=None)
        graph.add_pin(r'\SI{%s}{\degree}' % angle,
                      location=locations.next(), use_arrow=True,
                      relative_position=positions.next())

    xlabel("Core distance [m]")
    graph.set_xlabel(r"Core distance [\si{\meter}]")
    ylabel("Detection efficiency")
    graph.set_ylabel("Detection efficiency")
    #title(r"$N_{MIP} \geq %d$" % N)
    legend()
    graph.set_xlimits(0, 100)
    graph.set_ylimits(0, 1)

    utils.saveplot(N)
    artist.utils.save_graph(graph, suffix=N, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def boxplot_phi_reconstruction_results_for_MIP(group, N):
    table = group.E_1PeV.zenith_22_5

    figure()

    bin_edges = linspace(-180, 180, 18)
    x, r_dphi = [], []
    d25, d50, d75 = [], [], []
    for low, high in zip(bin_edges[:-1], bin_edges[1:]):
        rad_low = deg2rad(low)
        rad_high = deg2rad(high)
        query = '(min_n134 >= N) & (rad_low < reference_phi) & (reference_phi < rad_high)'
        sel = table.read_where(query)
        dphi = sel[:]['reconstructed_phi'] - sel[:]['reference_phi']
        dphi = (dphi + pi) % (2 * pi) - pi
        r_dphi.append(rad2deg(dphi))

        d25.append(scoreatpercentile(rad2deg(dphi), 25))
        d50.append(scoreatpercentile(rad2deg(dphi), 50))
        d75.append(scoreatpercentile(rad2deg(dphi), 75))
        x.append((low + high) / 2)

    fill_between(x, d25, d75, color='0.75')
    plot(x, d50, 'o-', color='black')

    xlabel(r"$\phi_{simulated}$ [deg]")
    ylabel(r"$\phi_{reconstructed} - \phi_{simulated}$ [deg]")
    #title(r"$N_{MIP} \geq %d, \quad \theta = 22.5^\circ$" % N)

    xticks(linspace(-180, 180, 9))
    axhline(0, color='black')
    ylim(-15, 15)

    utils.saveplot(N)

    graph = GraphArtist()
    graph.draw_horizontal_line(0, linestyle='gray')
    graph.shade_region(x, d25, d75)
    graph.plot(x, d50, linestyle=None)
    graph.set_xlabel(r"$\phi_\mathrm{sim}$ [\si{\degree}]")
    graph.set_ylabel(r"$\phi_\mathrm{rec} - \phi_\mathrm{sim}$ [\si{\degree}]")
    graph.set_title(r"$N_\mathrm{MIP} \geq %d$" % N)
    graph.set_xticks([-180, -90, '...', 180])
    graph.set_xlimits(-180, 180)
    graph.set_ylimits(-17, 17)
    artist.utils.save_graph(graph, suffix=N, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
    def plot_detection_efficiency(self):
        integrals, dens = self.get_integrals_and_densities()

        popt = self.full_fit_on_data(integrals,
                                      (1., 1., 5e3 / .32, 3.38 / 5000, 1.))

        x, y, yerr = [], [], []
        dens_bins = np.linspace(0, 10, 51)
        for low, high in zip(dens_bins[:-1], dens_bins[1:]):
            sel = integrals.compress((low <= dens) & (dens < high))
            x.append((low + high) / 2)
            frac = self.determine_charged_fraction(sel, popt)
            y.append(frac)
            yerr.append(np.sqrt(frac * len(sel)) / len(sel))
            print (low + high) / 2, len(sel)
            self.plot_full_spectrum_fit_in_density_range(sel, popt, low, high)
        print

        plt.figure()
        plt.errorbar(x, y, yerr, fmt='o', label='data', markersize=3.)

        popt, pcov = optimize.curve_fit(self.conv_p_detection, x, y, p0=(1.,))
        print "Sigma Gauss:", popt

        x2 = plt.linspace(0, 10, 101)
        plt.plot(x2, self.p_detection(x2), label='poisson')
        plt.plot(x2, self.conv_p_detection(x2, *popt), label='poisson/gauss')

        plt.xlabel("Charged particle density [$m^{-2}$]")
        plt.ylabel("Detection probability")
        plt.ylim(0, 1.)
        plt.legend(loc='best')
        utils.saveplot()

        graph = GraphArtist()
        graph.plot(x2, self.p_detection(x2), mark=None)
        graph.plot(x2, self.conv_p_detection(x2, *popt), mark=None,
                   linestyle='dashed')
        graph.plot(x, y, yerr=yerr, linestyle=None)
        graph.set_xlabel(
            r"Charged particle density [\si{\per\square\meter}]")
        graph.set_ylabel("Detection probability")
        graph.set_xlimits(min=0)
        graph.set_ylimits(min=0)
        artist.utils.save_graph(graph, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_front_passage():
    sim = data.root.showers.E_1PeV.zenith_0.shower_0
    leptons = sim.leptons
    R = 40
    dR = 2
    low = R - dR
    high = R + dR
    global t
    t = leptons.read_where('(low < core_distance) & (core_distance <= high)',
                          field='arrival_time')

    n, bins, patches = hist(t, bins=linspace(0, 30, 31), histtype='step')

    graph = GraphArtist()
    graph.histogram(n, bins)
    graph.set_xlabel(r"Arrival time [\si{\nano\second}]")
    graph.set_ylabel("Number of leptons")
    graph.set_ylimits(min=0)
    graph.set_xlimits(0, 30)
    graph.save('plots/front-passage')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def artistplot_reconstruction_efficiency_vs_R_for_angles(N):
    filename = 'DIR-plot_reconstruction_efficiency_vs_R_for_angles-%d.txt' % N
    all_data = loadtxt(os.path.join('plots/', filename))

    graph = GraphArtist()
    locations = iter(['above right', 'below left', 'below left'])
    positions = iter([.9, .2, .2])

    x = all_data[:, 0]

    for angle, efficiencies in zip([0, 22.5, 35], all_data[:, 1:].T):
        graph.plot(x, efficiencies, mark=None)
        graph.add_pin(r'\SI{%s}{\degree}' % angle, use_arrow=True,
                      location=locations.next(),
                      relative_position=positions.next())

    graph.set_xlabel("Core distance [\si{\meter}]")
    graph.set_ylabel("Reconstruction efficiency")
    graph.set_xlimits(0, 100)
    graph.set_ylimits(max=1)
    artist.utils.save_graph(graph, suffix=N, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
    def plot_gamma_landau_fit(self):
        events = self.data.root.hisparc.cluster_kascade.station_601.events
        ph0 = events.col('integrals')[:, 0]

        bins = np.linspace(0, RANGE_MAX, N_BINS + 1)
        n, bins = np.histogram(ph0, bins=bins)
        x = (bins[:-1] + bins[1:]) / 2

        p_gamma, p_landau = self.full_spectrum_fit(x, n, (1., 1.),
                                                   (5e3 / .32, 3.38 / 5000, 1.))
        print "FULL FIT"
        print p_gamma, p_landau

        n /= 10
        p_gamma, p_landau = self.constrained_full_spectrum_fit(x, n, p_gamma, p_landau)
        print "CONSTRAINED FIT"
        print p_gamma, p_landau

        plt.figure()
        print self.calc_charged_fraction(x, n, p_gamma, p_landau)

        plt.plot(x * VNS, n)
        self.plot_landau_and_gamma(x, p_gamma, p_landau)
        #plt.plot(x, n - self.gamma_func(x, *p_gamma))
        plt.xlabel("Pulse integral [V ns]")
        plt.ylabel("Count")
        plt.yscale('log')
        plt.xlim(0, 30)
        plt.ylim(1e1, 1e4)
        plt.legend()
        utils.saveplot()

        graph = GraphArtist('semilogy')
        graph.histogram(n, bins * VNS, linestyle='gray')
        self.artistplot_landau_and_gamma(graph, x, p_gamma, p_landau)
        graph.set_xlabel(r"Pulse integral [\si{\volt\nano\second}]")
        graph.set_ylabel("Count")
        graph.set_xlimits(0, 30)
        graph.set_ylimits(1e1, 1e4)
        artist.utils.save_graph(graph, dirname='plots')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
    def plot_full_spectrum_fit_in_density_range(self, sel, popt, low, high):
        bins = np.linspace(0, RANGE_MAX, N_BINS + 1)
        n, bins = np.histogram(sel, bins=bins)
        x = (bins[:-1] + bins[1:]) / 2

        p_gamma, p_landau = self.constrained_full_spectrum_fit(x, n, popt[:2], popt[2:])

        plt.figure()
        plt.plot(x * VNS, n, label="data")
        self.plot_landau_and_gamma(x, p_gamma, p_landau)

        y_charged = self.calc_charged_spectrum(x, n, p_gamma, p_landau)
        plt.plot(x * VNS, y_charged, label="charged particles")

        plt.yscale("log")
        plt.xlim(0, 50)
        plt.ylim(ymin=1)
        plt.xlabel("Pulse integral [V ns]")
        plt.ylabel("Count")
        plt.legend()
        suffix = "%.1f-%.1f" % (low, high)
        suffix = suffix.replace(".", "_")
        utils.saveplot(suffix)

        n = np.where(n > 0, n, 1e-99)
        y_charged = np.where(y_charged > 0, y_charged, 1e-99)

        graph = GraphArtist("semilogy")
        graph.histogram(n, bins * VNS, linestyle="gray")
        self.artistplot_alt_landau_and_gamma(graph, x, p_gamma, p_landau)
        graph.histogram(y_charged, bins * VNS)
        graph.set_xlabel(r"Pulse integral [\si{\volt\nano\second}]")
        graph.set_ylabel("Count")
        graph.set_title(
            r"$\SI{%.1f}{\per\square\meter} \leq \rho_\mathrm{charged}$ < $\SI{%.1f}{\per\square\meter}$" % (low, high)
        )
        graph.set_xlimits(0, 30)
        graph.set_ylimits(1e0, 1e4)
        artist.utils.save_graph(graph, suffix, dirname="plots")

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_R():
    graph = GraphArtist(width=r'.45\linewidth')

    n, bins, patches = hist(data.root.simulations.E_1PeV.zenith_22_5.shower_0.coincidences.col('r'), bins=100, histtype='step')
    graph.histogram(n, bins, linestyle='black!50')

    shower = data.root.simulations.E_1PeV.zenith_22_5.shower_0
    ids = shower.observables.get_where_list('(n1 >= 1) & (n3 >= 1) & (n4 >= 1)')
    R = shower.coincidences.read_coordinates(ids, field='r')
    n, bins, patches = hist(R, bins=100, histtype='step')
    graph.histogram(n, bins)

    xlabel("Core distance [m]")
    ylabel("Number of events")

    print("mean", mean(R))
    print("median", median(R))

    graph.set_xlabel(r"Core distance [\si{\meter}]")
    graph.set_ylabel("Number of events")
    graph.set_xlimits(min=0)
    graph.set_ylimits(min=0)
    graph.save('plots/SIM-R')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def boxplot_core_distance_vs_time():
    plt.figure()

    sim = data.root.showers.E_1PeV.zenith_0.shower_0
    leptons = sim.leptons

    #bins = np.logspace(0, 2, 25)
    bins = np.linspace(0, 100, 15)
    x, arrival_time, widths = [], [], []
    t25, t50, t75 = [], [], []
    for low, high in zip(bins[:-1], bins[1:]):
        sel = leptons.read_where('(low < core_distance) & (core_distance <= high)')
        x.append(np.mean([low, high]))
        arrival_time.append(sel[:]['arrival_time'])
        widths.append((high - low) / 2)
        ts = sel[:]['arrival_time']
        t25.append(scoreatpercentile(ts, 25))
        t50.append(scoreatpercentile(ts, 50))
        t75.append(scoreatpercentile(ts, 75))

    fill_between(x, t25, t75, color='0.75')
    plot(x, t50, 'o-', color='black')

    plt.xlabel("Core distance [m]")
    plt.ylabel("Arrival time [ns]")

    #utils.title("Shower front timing structure")
    utils.saveplot()

    graph = GraphArtist()
    graph.plot(x, t50, linestyle=None)
    graph.shade_region(x, t25, t75)
    graph.set_xlabel(r"Core distance [\si{\meter}]")
    graph.set_ylabel(r"Arrival time [\si{\nano\second}]")
    graph.set_ylimits(0, 30)
    graph.set_xlimits(0, 100)
    graph.save('plots/front-passage-vs-R')

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_pulseheight_histogram(data):
    events = data.root.hisparc.cluster_kascade.station_601.events
    ph = events.col("pulseheights")

    s = landau.Scintillator()
    mev_scale = 3.38 / 340
    count_scale = 6e3 / 0.32

    clf()
    n, bins, patches = hist(ph[:, 0], bins=arange(0, 1501, 10), histtype="step")
    x = linspace(0, 1500, 1500)
    plot(x, s.conv_landau_for_x(x, mev_scale=mev_scale, count_scale=count_scale))
    plot(x, count_scale * s.landau_pdf(x * mev_scale))
    ylim(ymax=25000)
    xlim(xmax=1500)

    # Remove one statistical fluctuation from data.  It is not important
    # for the graph, but it detracts from the main message
    index = bins.searchsorted(370)
    n[index] = mean([n[index - 1], n[index + 1]])

    graph = GraphArtist()
    n_trunc = where(n <= 100000, n, 100000)
    graph.histogram(n_trunc, bins, linestyle="gray")
    graph.add_pin("data", x=800, location="above right", use_arrow=True)
    graph.add_pin("$\gamma$", x=90, location="above right", use_arrow=True)
    graph.plot(x, s.conv_landau_for_x(x, mev_scale=mev_scale, count_scale=count_scale), mark=None)
    graph.add_pin("convolved Landau", x=450, location="above right", use_arrow=True)
    graph.plot(x, count_scale * s.landau_pdf(x * mev_scale), mark=None, linestyle="black")
    graph.add_pin("Landau", x=380, location="above right", use_arrow=True)

    graph.set_xlabel(r"Pulseheight [\adc{}]")
    graph.set_ylabel(r"Number of events")
    graph.set_xlimits(0, 1400)
    graph.set_ylimits(0, 21000)
    graph.save("plots/plot_pulseheight_histogram")

# 需要导入模块: from artist import GraphArtist [as 别名]
# 或者: from artist.GraphArtist import set_xlimits [as 别名]
def plot_uncertainty_mip(table):
    rec = DirectionReconstruction

    # constants for uncertainty estimation
    station = table.attrs.cluster.stations[0]
    r1, phi1 = station.calc_r_and_phi_for_detectors(1, 3)
    r2, phi2 = station.calc_r_and_phi_for_detectors(1, 4)

    THETA = deg2rad(22.5)
    DTHETA = deg2rad(5.)
    DN = .1
    LOGENERGY = 15
    DLOGENERGY = .5

    figure()
    x, y, y2 = [], [], []
    for N in range(1, 6):
        x.append(N)
        events = table.read_where('(abs(min_n134 - N) <= DN) & (abs(reference_theta - THETA) <= DTHETA) & (abs(log10(k_energy) - LOGENERGY) <= DLOGENERGY)')
        print(len(events),)
        errors = events['reference_theta'] - events['reconstructed_theta']
        # Make sure -pi < errors < pi
        errors = (errors + pi) % (2 * pi) - pi
        errors2 = events['reference_phi'] - events['reconstructed_phi']
        # Make sure -pi < errors2 < pi
        errors2 = (errors2 + pi) % (2 * pi) - pi
        #y.append(std(errors))
        #y2.append(std(errors2))
        y.append((scoreatpercentile(errors, 83) - scoreatpercentile(errors, 17)) / 2)
        y2.append((scoreatpercentile(errors2, 83) - scoreatpercentile(errors2, 17)) / 2)

    print()
    print("mip: min_n134, theta_std, phi_std")
    for u, v, w in zip(x, y, y2):
        print(u, v, w)
    print()

    # Simulation data
    sx, sy, sy2 = loadtxt(os.path.join(DATADIR, 'DIR-plot_uncertainty_mip.txt'))

    # Uncertainty estimate
    ex = linspace(1, 5, 50)
    phis = linspace(-pi, pi, 50)
    phi_errsq = mean(rec.rel_phi_errorsq(pi / 8, phis, phi1, phi2, r1, r2))
    theta_errsq = mean(rec.rel_theta1_errorsq(pi / 8, phis, phi1, phi2, r1, r2))
    #ey = TIMING_ERROR * std_t(ex) * sqrt(phi_errsq)
    #ey2 = TIMING_ERROR * std_t(ex) * sqrt(theta_errsq)

    R_list = [30, 20, 16, 14, 12]
    with tables.open_file('master-ch4v2.h5') as data2:
        mc = my_std_t_for_R(data2, x, R_list)
    mc = sqrt(mc ** 2 + 1.2 ** 2 + 2.5 ** 2)
    print(mc)
    ey = mc * sqrt(phi_errsq)
    ey2 = mc * sqrt(theta_errsq)

    nx = linspace(1, 4, 100)
    ey = spline(x, ey, nx)
    ey2 = spline(x, ey2, nx)

    # Plots
    plot(x, rad2deg(y), '^', label="Theta")
    plot(sx, rad2deg(sy), '^', label="Theta (sim)")
    plot(nx, rad2deg(ey2))#, label="Estimate Theta")
    plot(x, rad2deg(y2), 'v', label="Phi")
    plot(sx, rad2deg(sy2), 'v', label="Phi (sim)")
    plot(nx, rad2deg(ey))#, label="Estimate Phi")

    # Labels etc.
    xlabel("$N_{MIP} \pm %.1f$" % DN)
    ylabel("Angle reconstruction uncertainty [deg]")
    title(r"$\theta = 22.5^\circ \pm %d^\circ \quad %.1f \leq \log(E) \leq %.1f$" % (rad2deg(DTHETA), LOGENERGY - DLOGENERGY, LOGENERGY + DLOGENERGY))
    legend(numpoints=1)
    xlim(0.5, 4.5)
    utils.saveplot()
    print

    graph = GraphArtist()
    graph.plot(x, rad2deg(y), mark='o', linestyle=None)
    graph.plot(sx, rad2deg(sy), mark='square', linestyle=None)
    graph.plot(nx, rad2deg(ey2), mark=None)
    graph.plot(x, rad2deg(y2), mark='*', linestyle=None)
    graph.plot(sx, rad2deg(sy2), mark='square*', linestyle=None)
    graph.plot(nx, rad2deg(ey), mark=None)
    graph.set_xlabel(r"$N_\mathrm{MIP} \pm %.1f$" % DN)
    graph.set_ylabel(r"Angle reconstruction uncertainty [\si{\degree}]")
    graph.set_xlimits(max=4.5)
    graph.set_ylimits(0, 40)
    graph.set_xticks(range(5))
    artist.utils.save_graph(graph, dirname='plots')