How to use Landau - 3 common examples

from __future__ import division

import tables
from scipy import optimize

from landau import Scintillator

try:
    data
except NameError:
    data = tables.openFile('kascade.h5', 'r')

events = data.root.efficiency.events.read()
dens = events['k_cosdens_charged'][:,1]
ph0 = events[:]['pulseheights'][:,1]
s = Scintillator()

figure()
# Fit of convoluted Landau
n, bins, patches = hist(ph0, bins=linspace(0, 2000, 101), histtype='step',
                        label="Data")
nx = bins[:-1] + .5 * (bins[1] - bins[0])
x = linspace(-2000, 2000, 200)
y = interp(x, nx, n)
p = optimize.fmin(s.residuals, (10 ** 4, 3.38 / 380., 1), (x, y, 350, 500))
plot(x, s.conv_landau(x, *p), label='Charged particles')

# Fit of gamma spectrum
f = lambda x, N, a: N * x ** -a
x2 = x.compress((0 <= x) & (x < 100))
y2 = y.compress((0 <= x) & (x < 100))
popt, pcov = optimize.curve_fit(f, x2, y2, sigma=y2)

def plot_charged_particles_poisson(use_known=False):
    events = data.getNode(GROUP, 'events')
    s = Scintillator()

    if use_known:
        s.mev_scale = 0.0086306040898338834
        s.gauss_scale = 0.84289265239940525
        s.mev_scale *= 0.9
    else:
        ph = events[:]['pulseheights'][:,1]
        analyze_charged_particle_spectrum(s, ph, constrained=False)

    bins = linspace(0, 5, 41)
    x = bins[:-1] + .5 * (bins[1] - bins[0])
    y, yerr = [], []

    events = events.read()
    dens = events['k_cosdens_charged'][:,1]     # center detector
    for num, (u, v) in enumerate(zip(bins[:-1], bins[1:])):

figure()

    x = data.root.datasets.poisson.x.read()
    y = data.root.datasets.poisson.y.read()
    yerr = data.root.datasets.poisson.yerr.read()
    errorbar(x, y, yerr=yerr, fmt='o', label="Data")

    x = linspace(-10, 10, 101)
    f = vectorize(lambda x: 1 - exp(-.5 * x) if x >= 0 else 0.)

    for s in [0., 0.5, 1., 1.5]:
        if s == 0.:
            plot(x, f(x), label="Unconvoluted")
        else:
            g = stats.norm(scale=s).pdf
            plot(x, discrete_convolution(f, g, x),
                 label="sigma = %.2f m$^{-2}$" % s)

    xlim(0, 5)
    title("Effect of uncertainty on particle density")
    xlabel("Charged particle density (m$^{-2}$)")
    ylabel("Probability of one or more particles")
    legend(loc='best')
    savefig("plots/conv_poisson.pdf")