fit

fit.py

#!/usr/bin/env python2.7

import numpy
import scipy
import scipy.optimize
import subprocess


def plotme(gnuplottext, filename=None):
    if filename:
        fout = open(filename, "w")
        fout.write(gnuplottext)
        fout.close()
        gnusubprocess = subprocess.Popen(['gnuplot', filename])
    else:
        gnusubprocess = subprocess.Popen(['gnuplot'], STDIN=gnuplottext)


def gauss_func(p, x, y):
    fx = p[0] * numpy.exp(-(x-p[1])**2/2.0/p[2]**2)+p[3] * numpy.exp(-(x-p[4])**2/2.0/p[5]**2)+p[6]
    residual = y - fx
    return residual

for filename in ('01.mca', '02.mca', '03.mca'):
    print "-----"
    print filename
    signal = []
    for line in open(filename):
        signal.append(int(line))

    inich1 = 168
    inich2 = 188
    inisig1 = 5.2
    inisig2 = 7.0
    ini_param = [signal[inich1], inich2, inisig1, signal[inich2], inich2, inisig2, 200]
    Px = numpy.array(range(int(inich1-8*inisig1), int(inich2+5*inisig2)))
    Py = numpy.array(signal[int(inich1-8*inisig1):int(inich2+5*inisig2)])

    result, error = scipy.optimize.leastsq(gauss_func, ini_param, args=(Px, Py))
    amp1, ch1, sig1, amp2, ch2, sig2, bg1 = result
    print "amp1: %12.4f" % amp1
    print "ch1:  %12.4f" % ch1
    print "sig1: %12.4f" % sig1
    print "amp2: %12.4f" % amp2
    print "ch2:  %12.4f" % ch2
    print "sig2: %12.4f" % sig2
    print "bg1:  %12.4f" % bg1

    gnutext = """
#!/usr/local/bin/gnuplot
set terminal aqua
set xrange[0:250]
set sample 1000
fitfunc(x) = (%12.4f)*exp(-(x-(%12.4f))**2/2.0/(%12.4f)**2)+(%12.4f)*exp(-(x-(%12.4f))**2/2.0/(%12.4f)**2)+(%12.4f)
""" % (amp1, ch1, sig1, amp2, ch2, sig2, bg1) + """
set title "%s"
""" % (filename) + """
plot "%s" with histeps
replot fitfunc(x)
""" % filename
    plotme(gnutext, filename=filename + ".plt")