Warning, /analysis/Prototype2/EMCal/ShowerCalib/Fit/Proto2ShowerCalibFit.m is written in an unsupported language. File is not indexed.
0001
0002 clear all
0003 close all
0004
0005 global Ndata;
0006 global DataSet;
0007 global low_tower_IDs;
0008
0009 %%
0010
0011 DataFolder = 'E:/tmp/ShowerCalib/';
0012
0013 % FileID = {'Rot45','THP','UIUC18','UpTilt5', 'ShowerDepth'};
0014 FileID = {'Rot45','THP','UIUC18','UIUC21','Tilt0', 'ShowerDepth'};
0015 FileList = FileID;
0016
0017 sim_const = 3/100;
0018 sim_stat = 12/100;
0019 Ndata = 64;
0020 InitConst = [ones(1, 64) ];
0021
0022 zero_sup = -1000; %zero suppression disabled
0023
0024 %%
0025
0026 for i = 1:size(FileList,2)
0027
0028 FileList{i} = [DataFolder FileList{i} '.lst_EMCalCalib.root.dat'];
0029
0030 end
0031
0032 %%
0033
0034 % N_Runs = size(RunList, 2);
0035
0036 %% Readin
0037 % global DataSet
0038 DataSet =struct('FileID',{},'E',{},'DE',{},'data',{}, 'accept',{});
0039
0040 N_Runs = 0;
0041 % for i = 1:1
0042 for i = 1:size(FileList,2)
0043 filename = FileList{i};
0044 fprintf('processing %s\n', filename);
0045
0046 data = textread(filename);
0047 % disp(size(data));
0048
0049 energys = data(:,1);
0050 energy = unique(energys);
0051 data = data(:,2:(1+Ndata));
0052
0053 energy = energy(energy<=8 & energy>2);
0054
0055 for j=1:size(energy, 1)
0056
0057 N_Runs = N_Runs +1;
0058 fprintf('processing %s - %.0f GeV @ %d\n', filename,energy(j), N_Runs );
0059
0060 DataSet(N_Runs).FileID = FileID{i};
0061 DataSet(N_Runs).E = energy(j);
0062 DataSet(N_Runs).DE = sqrt( sim_const.^2 + sim_stat.^2./DataSet(i).E ) ;
0063
0064 DataSet(N_Runs).data = data(energys == energy(j),:);
0065 DataSet(N_Runs).data = DataSet(N_Runs).data .* (DataSet(N_Runs).data>zero_sup);
0066
0067 total_E = sum( DataSet(N_Runs).data* InitConst', 2) ;
0068
0069 DataSet(N_Runs).data = DataSet(N_Runs).data(total_E > DataSet(N_Runs).E/2, :);
0070 DataSet(N_Runs).accept = ones(size(data, 1), 1);
0071
0072 end
0073
0074 end
0075
0076 %%
0077 DrawDataSet(DataSet, InitConst, 'Inputs');
0078
0079 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0080
0081
0082
0083 %%
0084
0085 figure('name',['DrawDataSet_EnergyFraction'],'PaperPositionMode','auto', ...
0086 'position',[100,0,1800,1000]) ;
0087
0088 SumFraction = [];
0089
0090 for i = 1:N_Runs
0091 subplot(4,4,i);
0092 total_E = ones(size(DataSet(i).data)) .*DataSet(i).E;
0093 Fraction = DataSet(i).data ./ total_E;
0094 MeanFraction = mean(Fraction, 1);
0095 MeanFraction = reshape(MeanFraction, 8, 8);
0096
0097 SumFraction = [SumFraction; Fraction];
0098
0099 imagesc(0:7, 0:7, MeanFraction);
0100 colorbar
0101 set(gca,'YDir','normal')
0102
0103 title(sprintf('%s, E = %.1f GeV', DataSet(i).FileID, DataSet(i).E));
0104 xlabel('Column ID');
0105 ylabel('Row ID');
0106 end
0107
0108 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0109 %%
0110
0111 figure('name',['DrawDataSet_EnergyFractionSum'],'PaperPositionMode','auto', ...
0112 'position',[100,0,1900,400]) ;
0113
0114 subplot(1,3,1);
0115
0116 MeanFraction = mean(SumFraction, 1);
0117 MeanFraction = reshape(MeanFraction, 8, 8);
0118
0119 imagesc(0:7, 0:7, MeanFraction);
0120 colorbar
0121 set(gca,'YDir','normal');
0122
0123 title(sprintf('Sum all data'));
0124 xlabel('Column ID');
0125 ylabel('Row ID');
0126
0127 subplot(1,3,2);
0128
0129 hist(reshape(MeanFraction,size(MeanFraction,1) * size(MeanFraction,2),1),100);
0130 title('fraction of shower energy carried by each tower');
0131 xlabel('MeanFraction');
0132
0133 subplot(1,3,3);
0134
0135 low_tower_IDs =reshape( MeanFraction<0.004, 1, Ndata);
0136
0137 imagesc(0:7, 0:7, reshape(low_tower_IDs, 8, 8));
0138 colorbar
0139 set(gca,'YDir','normal');
0140
0141 title(sprintf('Low hit towers'));
0142 xlabel('Column ID');
0143 ylabel('Row ID');
0144
0145 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0146
0147 % return
0148 %%
0149
0150 InitConst_RunScale = [InitConst ones(1, N_Runs)];
0151 % InitConst_RunScale = [InitConst ];
0152
0153 data_selection(InitConst_RunScale,10);
0154 disp(object_function(InitConst_RunScale));
0155
0156 % disp(object_function(InitConst_RunScale, 2));
0157 % disp(object_function(InitConst_RunScale, 1));
0158
0159 %%
0160
0161 x = InitConst_RunScale;
0162
0163 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0164
0165 data_selection(x,8);
0166 disp(object_function(x));
0167
0168 x = fminsearch(@(x) object_function(x), x,...
0169 options);
0170
0171 %%
0172 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0173 %
0174 data_selection(x,4);
0175 disp(object_function(x));
0176
0177 x = fminsearch(@(x) object_function(x), x,...
0178 options);
0179
0180
0181 %%
0182 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0183 %
0184 data_selection(x,2);
0185 disp(object_function(x));
0186
0187 x = fminsearch(@(x) object_function(x), x,...
0188 options);
0189
0190
0191 %%
0192
0193 data_selection(x,2);
0194
0195 calib_const = x(1:Ndata);
0196 E_scale = x((Ndata+1):(Ndata + N_Runs));
0197
0198
0199 figure('name',['CalibConst'],'PaperPositionMode','auto', ...
0200 'position',[100,0,1800,400]) ;
0201
0202 subplot(1,3,1);
0203
0204
0205 plot(calib_const);
0206 title(sprintf('Calibration constant'));
0207 xlabel('Col * 8 + Row');
0208 ylabel('Calibration New / Old');
0209
0210 subplot(1,3,2);
0211
0212 imagesc(0:7, 0:7, reshape(calib_const, 8, 8));
0213 colorbar
0214 set(gca,'YDir','normal')
0215
0216 title(sprintf('Calibration constant, New / Old'));
0217 xlabel('Column ID');
0218 ylabel('Row ID');
0219
0220
0221 subplot(1,3,3);
0222
0223 plot(E_scale);
0224 title(sprintf('Energy scale constant, mean = %.1f', mean(E_scale)));
0225 xlabel('Run ID');
0226 ylabel('Energy scale New / Old');
0227
0228
0229 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0230 %%
0231
0232 DrawDataSet(DataSet,calib_const,'Optimized');
0233 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0234
0235 %% Save out
0236 %
0237 % for i = 1:N_Runs
0238 % filename = sprintf('%s/calirbated_%d_%.0f.dat',DataFolder, DataSet(i).FileID, DataSet(i).E );
0239 %
0240 %
0241 % total_E = sum( DataSet(i).data*InitConst', 2) ;
0242 % calib_total_E = sum( DataSet(i).data* calib_const', 2) ;
0243 %
0244 % dlmwrite(filename,[total_E calib_total_E]);
0245 % end
0246
0247 save([DataFolder 'fit.mat']);
0248 % save('goodfit.mat');
0249
0250 %%
0251
0252 [calib_const_col, calib_const_row]= meshgrid(0:7,0:7);
0253
0254 A = [reshape(calib_const_col,1,64); reshape(calib_const_row,1,64); calib_const];
0255
0256 fileID = fopen([DataFolder 'ShowerCalibFit_CablibConst.dat'],'w');
0257 % fprintf(fileID,'%d - %d \n',);
0258
0259 fprintf(fileID,'%d\t%d\t%f\n',A);
0260
0261
0262 fclose(fileID);
0263
0264
0265
0266
