1. 程式人生 > >MATLAB之TOPAS輸出檔案擬合高斯分佈並提取引數

MATLAB之TOPAS輸出檔案擬合高斯分佈並提取引數

此處需要將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來跑。