EMCal/macros/EnergyCalibFit.m

0001
0002 clear all
0003 close all
0004
0005 %%
0006
0007 %   const double Es[] =
0008 %     { 2,   3  ,  4,    8,   12 ,   16 };
0009 %   const double runs[] =
0010 %     { 2042,2040, 2039, 2038, 2067, 2063 };
0011
0012 RunList = [
0013     2,   3  ,  4,    8,   12 ,   16
0014     2042,2040, 2039, 2038, 2067, 2063
0015 ];
0016 N_Runs = size(RunList, 2);
0017
0018 sim_const = 2.409/100;
0019 sim_stat = 11.77/100;
0020 Ndata = 20;
0021
0022 %% Readin
0023 % global DataSet
0024 DataSet =struct('run',{},'E',{},'DE',{},'data',{});
0025
0026 for i = 1:N_Runs
0027
0028     DataSet(i).run = RunList(2,i);
0029     DataSet(i).E = RunList(1,i);
0030     DataSet(i).DE = sqrt( sim_const.^2 + sim_stat.^2./DataSet(i).E  )  *  DataSet(i).E;
0031
0032     filename = sprintf('beam_0000%d-0000_DSTReader.root_DrawPrototype2EMCalTower_Prototype2_DSTReader.dat', DataSet(i).run );
0033
0034     fprintf('processing %s\n', filename);
0035
0036     data = textread(filename);
0037     data = data(:,1:Ndata);
0038
0039     total_E = sum(data, 2);
0040     data = data(total_E>(DataSet(i).E * 0.5) & total_E<(DataSet(i).E * 1.5) , :);
0041
0042     DataSet(i).data =  data;
0043
0044 end
0045
0046 %%
0047 DrawDataSet(DataSet, ones(1, Ndata), 'Inputs');
0048
0049 SaveCanvas(['EnergyCalibFIt'],gcf);
0050 %%
0051
0052 disp(object_function(ones(1, Ndata + N_Runs -1),DataSet, 10));
0053 disp(object_function(ones(1, Ndata + N_Runs -1),DataSet, 2));
0054 disp(object_function(ones(1, Ndata + N_Runs -1),DataSet, 1));
0055
0056 %%
0057
0058 options = optimset('Display','iter','TolFun',10000, 'MaxFunEvals', 100000,'MaxIter',100000,'PlotFcns',@optimplotfval );
0059
0060 x = fminsearch(@(x) object_function(x,DataSet, 10),ones(1, Ndata + N_Runs -1),...
0061     options);
0062
0063 %%
0064 options = optimset('Display','iter','TolFun',1, 'MaxFunEvals', 100000,'MaxIter',100000,'PlotFcns',@optimplotfval );
0065
0066 x = fminsearch(@(x) object_function(x,DataSet,3),x,...
0067     options);
0068
0069 %%
0070 options = optimset('Display','iter','TolFun',1e-4, 'MaxFunEvals', 100000,'MaxIter',100000,'PlotFcns',@optimplotfval );
0071
0072 x = fminsearch(@(x) object_function(x,DataSet,2),x,...
0073     options);
0074
0075 %%
0076
0077 calib_const = abs(x(1:Ndata));
0078
0079 DrawDataSet(DataSet,calib_const,'Optimized');
0080 SaveCanvas(['EnergyCalibFIt'],gcf);
0081
0082 %% Save out
0083
0084 for i = 1:N_Runs
0085     filename = sprintf('calirbated_%d.dat', DataSet(i).run );
0086
0087
0088     total_E = sum( DataSet(i).data, 2) ;
0089     calib_total_E = sum( DataSet(i).data* calib_const', 2) ;
0090
0091     dlmwrite(filename,[total_E calib_total_E]);
0092 end
0093
0094 save('fit.mat');
0095 % save('goodfit.mat');