此处需要将TOPAS跑出来的数据进行高斯拟合,并计算标准差sigma,以准备下一步的分析使用。csv数据类型如上篇文章所述,因为对TOPAS程序直接设定了输出文件的格式是“能量MeV+束斑大小cm+RangeShift宽度x空气间隙宽度”,此处用循环读取程序并输出到一个保存结果的xlsx文件中。本来可以直接用list=dir()函数以及list(i).name来跑程序,可是给我的文件命名格式有点问题,所以得重写循环。MATLAB用dir读取文件名的时候给出的排序都是从第一位开始依次比较大小,所以一开始定好文件名格式对后面的程序帮助很大。
clc;
clear;
fileDir='C:\Users\32628\Documents\MATLAB\results\';
%创建用于存放Sigma的数组
valueM=zeros(10,18);
row=1;
column=1;
for k=70:10:100
for j=1:1:3+(k-70)/10
for i=5:5:50
%读取文件名
kk=num2str(k);
jj=num2str(j);
ii=num2str(i);
filename=strcat(kk,'MeV',jj,'cm',ii,'.csv');
filename2=strcat(fileDir,filename);
file=csvread(filename2,8,0);
yvalue=file(:,4);
xvalue=file(:,1);
%设置要用的拟合模型
[xData, yData] = prepareCurveData( xvalue, yvalue );
ft = fittype( 'gauss1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
%opts.Display = 'Off';
%opts.Lower = [-Inf -Inf 0];
%opts.StartPoint = [39194592 299 24.6030834338969];
%将模型用于数据
[fitresult, gof] = fit( xData, yData, ft, opts );
% Plot fit with data.
%figure( 'Name', 'Normal Fit' );
%h = plot( fitresult, xData, yData );
%legend( h, 'yvalue vs. xvalue', 'untitled fit 1', 'Location', 'NorthEast' );
% Label axes
%xlabel xvalue
%ylabel yvalue
%grid on
valueM(row,column)=(fitresult.c1)/2^(1/2);
if i~=50
row=row+1;
else
row=1;
column=column+1;
end
end
end
end
xlswrite('C:\Users\32628\Documents\MATLAB\sigma.xlsx',valueM,'Sheet1','B2');
本来想直接用normfit()来得到sigma的,但发现用这个函数得到的sigma有问题,才另想办法用cftool来拟合数据并用Generate Code来跑。