File: lorentzian.py

package info (click to toggle)
python-bumps 1.0.3-3
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 6,200 kB
  • sloc: python: 24,517; xml: 493; ansic: 373; makefile: 211; javascript: 99; sh: 94
file content (105 lines) | stat: -rwxr-xr-x 2,871 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#!/usr/bin/env python

"""
Lorentzian peak model

References::
    None
"""

from numpy import array, pi, asarray, arange, sqrt
from numpy import random


def lorentzian(E, Io, Eo, Gamma, A=0, B=0, C=0):
    """
    lorentzian with quadratic background::

        I = Io  (Gamma/2*pi) / ( (E-Eo)^2 + (Gamma/2)^2 ) + (A + B E + C E^2)
    """
    E = asarray(E)  # force E to be a numpy array
    return A + (B + C * E) * E + Io * (Gamma / 2 / pi) / ((E - Eo) ** 2 + (Gamma / 2) ** 2)


def simulate_events(params, xmin, xmax, npts=4000):
    """Generate a lorentzian dataset of npts between [min,max] from given params"""

    def gensample(F, xmin, xmax):
        a = arange(xmin, xmax, (xmax - xmin) / 200.0)
        ymin = 0
        ymax = F(a).max()
        while 1:
            t1 = random.random() * (xmax - xmin) + xmin
            t2 = random.random() * (ymax - ymin) + ymin
            t3 = F(t1)
            if t2 < t3:
                return t1

    fwd = lambda x: lorentzian(x, **params)
    return array([gensample(fwd, xmin, xmax) for _ in range(npts)])


def simulate_histogram(pars, Emin, Emax, dE, npts=4000):
    events = simulate_events(pars, Emin, Emax)
    E, I = histogram(events, dE, Emin, Emax)
    # print min(events),max(events)
    dI = sqrt(I)
    data = {"x": (E[1:] + E[:-1]) / 2.0, "y": I, "dy": dI}
    return data


def demo_data():
    # integrated intensity = 4000
    # center = 6400
    # width gamma = 180
    # background = -(E-6340)^2/1000 + 10
    bgC = 6340.0
    bgW = 10000.0
    A = 35 - bgC**2 / bgW
    B = 2.0 * bgC / bgW
    C = -1.0 / bgW
    # A,B,C = 0,0,0
    N = 4000
    pars = {"Eo": 6500.0, "Gamma": 180.0, "Io": 20 * N, "A": A, "B": B, "C": C}
    Emin, Emax = 6000, 6700
    dE = 20

    return pars, simulate_histogram(pars, Emin, Emax, dE, N)


# probably shouldn't be in here...
from numpy import histogram as numpyhisto


def histogram(data, binwidth, xmin, xmax):
    """
    generate bin-centered histogram of provided data

    return bins of given binwidth (and histogram) generated between [xmin,xmax]
    """
    edges = arange(xmin, xmax + binwidth * 0.9999999, binwidth)
    centers = edges + (0.5 * binwidth)
    histo, _ = numpyhisto(data, bins=edges)
    # print data.size, sum(histo), edges[0], edges[-1], min(data),max(data)
    return centers, histo


coeff, data = demo_data()


def demo():
    import pylab

    x, y, dy = data["x"], data["y"], data["dy"]
    A, B, C = coeff["A"], coeff["B"], coeff["C"]
    Io, Eo, Gamma = coeff["Io"], coeff["Eo"], coeff["Gamma"]
    pylab.errorbar(x, y, yerr=dy, label="data")
    pylab.plot(x, pylab.polyval([C, B, A], x), label="background")
    pylab.plot(x, lorentzian(x, Eo=Eo, Io=Io, Gamma=Gamma), label="peak")
    pylab.plot(x, lorentzian(x, **coeff), label="peak+bkg")
    pylab.legend()
    pylab.show()


if __name__ == "__main__":
    demo()