In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [300]:

    

range_data = np.array([0.8516, 1.542, 2.866, 5.70, 9.15, 26.76, 36.96, 58.79, 93.32,
                152.4, 211.5, 441.8, 553.4, 771.2, 1088, 1599, 2095, 3998,
                4920, 6724, 9360, 13620, 17760, 33430, 40840, 54950, 74590,
                104000, 130200, 212900,
                2.453e5,
                2.990e5,
                3.616e5,
                4.384e5,
                4.957e5,
                6.400e5,
                6.877e5,
                7.603e5,
                8.379e5,
                9.264e5,
                9.894e5,
                1.141e6,
                1.189e6])

momentum_data = np.array([47.04, 56.16, 68.02, 85.1, 100, 152.7, 176.4, 221.8,
                    286.8, 391.7, 494.5, 899.5, 1101, 1502, 2103, 3104, 4104,
                    8105, 10110, 14110, 20110, 30110, 40110, 80110, 100100,
                    140100, 200100, 300100, 400100, 800100,
                    1e6,
                    1.4e6,
                    2e6,
                    3e6,
                    4e6,
                    8e6,
                    1e7,
                    1.4e7,
                    2e7,
                    3e7,
                    4e7,
                    8e7,
                    1e8])


    

In [305]:

    

plot(np.log(range_data), np.log(momentum_data), 'k.')
xlim((-2, 15))
xlabel("log range")
ylabel("log momentum")


    

    
Out[305]:





<matplotlib.text.Text at 0x1430ac10>






    

In [222]:

    

# find the slope and intercept for the low range data


    

In [345]:

    

slope = (log(momentum_data[3]) - log(momentum_data[0])) / (log(range_data[3]) - log(range_data[0]))
slope


    

    
Out[345]:





0.31183384474750259

In [346]:

    

rint = log(momentum_data[0]) - slope * log(range_data[0])
rint


    

    
Out[346]:





3.9010907766258351

In [308]:

    

log(momentum_data[1])


    

    
Out[308]:





4.0282047597175552

In [338]:

    

rng = np.linspace(np.log(range_data[0]), np.log(range_data[-1]), 2000)


    

In [323]:

    

#def fn(x, a):
#    return np.sqrt(a * (x - 3.8509983035903717)) - 0.16063834596281262
def fn(x, a, b, c):
    res = a * (x + b) ** 2 + c
    #res = sqrt(a * abs(x - b)) + c
    return res


    

In [327]:

    

plot(np.log(range_data), np.log(momentum_data), 'k.')
ylim((-2, 15))
from scipy.optimize import curve_fit
popt, _ = curve_fit(fn, np.log(range_data), np.log(momentum_data))
plot(rng, fn(rng, *popt))


    

    
Out[327]:





[<matplotlib.lines.Line2D at 0x15c9af10>]






    

In [325]:

    

lres = fn(np.log(range_data), *popt) - np.log(momentum_data)
plot(np.log(momentum_data), lres, 'k.')


    

    
Out[325]:





[<matplotlib.lines.Line2D at 0x157b7a90>]






    

In [333]:

    

range_data[3]


    

    
Out[333]:





5.7000000000000002

In [348]:

    

a, b, c = popt
a, b, c


    

    
Out[348]:





(0.066278146047421563, 0.01506720229624511, 4.2081356458976158)

In [335]:

    

def newfit(x):
    a, b, c = popt
    res = a * (x + b) ** 2 + c
    low_idx = x < log(range_data[3])
    res[low_idx] = slope * x[low_idx] + rint
    return res


    

In [339]:

    

loglog(range_data, momentum_data, 'k.')
#ylim((-2, 4))
#xlim((3,6))
plot(np.exp(rng), np.exp(newfit(rng)))
xlabel("Range (g/cm^2)")
ylabel("Momentum (MeV/c)")
title("Fit to Groom 2001 Muon Range data")


    

    
Out[339]:





<matplotlib.text.Text at 0x1692f290>






    

In [ ]: