python ガウス関数の積分 区分求積

python

import pylab as py import numpy as np from scipy import optimize from scipy.stats import stats from pylab import plot import matplotlib.pyplot as plt import pandas as pd import matplotlib.lines as mlines from scipy.integrate import quad FrontPFT='Minor.csv'#csv file renamed so that the plot can be identoified for which CR-39 it has drawn#in the plot T091 will appear as ID df = pd.read_csv(FrontPFT, skip_blank_lines=True ,skiprows=1,header=None,skipfooter=0,engine='python') minor = df.iloc[:,0] draw=minor#[r>=0.5]#binwidth=0.48bi = np.linspace(0, 48, 300)data = plt.hist(draw.dropna(), bins = bi) m = 32.78s = 0.47 #first Gaussian function for the first peak (peaks counted from the right)def f(x, a, b, c):# return (a * py.exp(-(x - b)**2.0 / (2 *c**2))) #a= m = 32.78 s = 0.47 a = 3202 return (a * py.exp(-(x - m)**2.0 / (2 *s**2))) # first function fitting, x and y values from the array minorx = [0.5 * (data[1][i] + data[1][i+1]) for i in range(len(data[1])-1)]y = data[0] #cryve fitting instructionspopt, pcov = optimize.curve_fit(f, x, y)chi, p = stats.chisquare(popt)#this line calculate chi_s/ndf when m and s is automaticchi_square = (np.sum(((draw - m)**2)/(m)))*1/(len(draw)-2)#when m and s defined manually then chi_s/ndf should calculate this wayx_fit = py.linspace(m-5, m+5, 80000)y_fit = f(x_fit, *popt)plot(x_fit, y_fit, lw=3, color="r",ls="--")#this line will draw a fitted red line of a peak specified by b and c value #integration of first function f#gaussF = lambda x, a: f(x, a, m, s)f1 = lambda x : f(x, 10187, m, s)#interagting the function f for a=10187 for carbon peak, mean and sigma valuesresult = quad(f1,(m-3*s),(m+3*s))#number of entities within this limit in fitted curvenumPar = result[0]#1st entity from result, totla number of particlern= minor[(minor>(m-3*s)) & (minor<(m+3*s))] #this line shows the number of particle within +-3s in csv file print ("Numbe of particel based on real data= ",rn.shape )bs = mlines.Line2D([], [], label='ID: =%s \n $\mu$=%.1f \n $\sigma$=%.1f \n $\chi^2$/ndf = %.1f\n Total etch pits(fitted red line) = %.0f' %(str(FrontPFT[0:4]),m, s,chi_square,numPar))plt.legend(handles = [bs], loc="upper left", fontsize=12)plt.rcParams["figure.figsize"] = [9,9]plt.xlabel("Minor axis (um)", color="m", size=16)plt.ylabel("Count/bin", color="m", size=20)plt.xlim(5,40)#plt.ylim(0,300)#plt.title("Energy (MeV/n)=%0.2f, LET (MeV/mm)=%0.2f, Number =%0.0f" %(E,LET,N),size=24)plt.tick_params(axis='both', which='major', labelsize=16)manager = plt.get_current_fig_manager()#plt.show()manager.window.showMaximized()fig1=plt.savefig("_Avg.jpg")