clc; clear me = 0.511; % electron rest mass in MeV %En = 14; %enter beam energy here in MeV unit. En = 12.74; %enter beam energy here in MeV unit. %I_max_neg = -10; %enter smallest negative scan current. Here -7 Amps %I_max_pos = 7; %enter biggest positive scan current. Here 10 Amps I_start=-8; % Scan starting current %I_end=-2.0; % Scan ending current I_increment=0.5; scan_step=13; % do 41 scan. for ii = 1:scan_step I(ii) = I_start + (ii-1)*I_increment; end p = sqrt(En*En - me*me)/1000;% e- momentum in GeV/c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %enter the distance from the center of the quad to the screen S12=3.1;%(in m unit) S12(q4)=1.875 cm. S12(q1)=3.10 cm erS12 = 0.005; % assume 51 mm error. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %enter calibration here: calibration_x=0.04308*0.001; %m/px calibration_y=0.04228*0.001; %m/px %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % now we do the action...................................... disp('Read image ...') % get image %rehash path rmpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q4_Scan\42mA\positive\2_Amp'); %addpath(genpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q4_Scan\42mA\negatvie\2.4')); addpath('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q4_Scan\42mA\positive\2_Amp'); %cd('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q4_Scan\42mA\positive\2_Amp'); bg1 = imread('bg\1.bmp'); bg2 = imread('bg\2.bmp'); bg3 = imread('bg\3.bmp'); bg4 = imread('bg\4.bmp'); bg5 = imread('bg\5.bmp'); bg6 = imread('bg\6.bmp'); bg7 = imread('bg\7.bmp'); bg8 = imread('bg\8.bmp'); bg9 = imread('bg\9.bmp'); bg10 = imread('bg\10.bmp'); bg = imread('bg\1.bmp'); %bg=(bg1+bg2+bg3+bg4+bg5+bg6+bg7+bg8+bg9+bg10)/10; imd_bg = double(bg(:,:,1)); %im = imread('C:\Documents and Settings\setisadi\My Documents\MATLAB\analysis\Mar_18_2011\Q1_Scan\42mA\Matlab_Files\1.jpg'); % reading 1.png sg = imread('2.bmp'); % reading 1.png imd_sg = double(sg(:,:,1)); pic = (imd_sg-imd_bg); %imd = imd_sg; %imd = imd_bg; % determine image size in pixels dim=size(pic); %dim=size(imd); dim_x = dim(1); dim_y = dim(2); xunit='m'; yunit='m'; % profiles profx = sum(pic,1);%how to make projection profy = sum(pic,2)';%how to make projection xbins=size(profx); ybins=size(profy); x=(1:1:xbins(2)); y=(1:1:ybins(2)); % gauss fit with narrowed range %xl = 150; %xr = 550; %yl = 100; %yr = 480; xl = 1; xr = 766; yl = 1; yr = 576; %xl = 100; %xr = 650; %yl = 1; %yr = 576; X = x; Y = profx; %X = y; %Y = profy; plot(X',Y','b') hold on; % a(1) = base % a(2) = A; amplitude % a(3) = x0; center % a(4) = D; sigma_0 % a(5) = N; a1 = 100; a2 = 7000; a3 = 280; a4 = 1.25; a5 = 1.5; %SupGau = a(1)+a(2)*exp(-abs(X-a(3))/(a(4))).^a(5); % super gaussian a0 = [a1,a2,a3,a4,a5]; opts = optimset('TolX',1e-4,'MaxFunEvals',10000,'MaxIter',10000,'Display','on'); [fitpara,fval,flag,err] = fminsearch(@SupGau_devsum,a0,opts,X,Y); base = fitpara(1) amplitude = fitpara(2) peak_center = fitpara(3) sig0 = fitpara(4) N = fitpara(5) sig = sig0*(pi/2).^(2/fitpara(5)-1) %y_fit = 1/(sqrt(2*pi)*fitpara(1))*exp(-abs(X-fitpara(3)).^fitpara(2)/(2*fitpara(1).^fitpara(2))); y_fit = fitpara(1)+fitpara(2)*exp( -0.5*(abs(X-fitpara(3))/(fitpara(4))).^fitpara(5) ); plot(X',y_fit','r'); %The function fminsearch can be tuned by additional options, one can demand a certain accuracy, restrict the number of optimization steps etc. %For a detailed description consult Help