Warning, /analysis/Prototype4/EMCal/ShowerCalib/Fit/Proto4ShowerCalibFit.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/Transfer Buffer/ShowerCalib/';
0012
0013 % FileID = {'Rot45','THP','UIUC18','UpTilt5', 'ShowerDepth'};
0014 % FileID = {'dst'};
0015 % FileID = {'ShowerCalib','ShowerCalib_tilted'};
0016
0017
0018 FileID = {'5thPositionScan_dst','6thPositionScan_dst'};
0019 FileList = FileID;
0020
0021 energy_scale = 0.55*8/3000;
0022 fixed_energy = 8;
0023
0024
0025 sim_const = 3/100;
0026 sim_stat = 12/100;
0027 Ndata = 64;
0028 InitConst = [ones(1, 64) ];
0029
0030 zero_sup = -1000; %zero suppression disabled
0031
0032 %%
0033
0034 for i = 1:size(FileList,2)
0035
0036 FileList{i} = [DataFolder FileList{i} '.lst_EMCalCalib.root.dat'];
0037
0038 end
0039
0040 %%
0041
0042 % N_Runs = size(RunList, 2);
0043
0044 %% Readin
0045 % global DataSet
0046 DataSet =struct('FileID',{},'E',{},'DE',{},'data',{}, 'accept',{});
0047
0048 N_Runs = 0;
0049 % for i = 1:1
0050 for i = 1:size(FileList,2)
0051 filename = FileList{i};
0052 fprintf('processing %s\n', filename);
0053
0054 data = textread(filename);
0055 % disp(size(data));
0056
0057 data(:,1) = ones(size(data(:,1))) * fixed_energy; % force fix beam energy settings.
0058
0059 energys = data(:,1);
0060 energy = unique(energys);
0061 data = data(:,2:(1+Ndata));
0062
0063 data = data * energy_scale;
0064
0065 energy = energy(energy<=8 & energy>2);
0066
0067 for j=1:size(energy, 1)
0068
0069 N_Runs = N_Runs +1;
0070 fprintf('processing %s - %.0f GeV @ %d\n', filename,energy(j), N_Runs );
0071
0072 DataSet(N_Runs).FileID = FileID{i};
0073 DataSet(N_Runs).E = energy(j);
0074 DataSet(N_Runs).DE = sqrt( sim_const.^2 + sim_stat.^2./DataSet(i).E ) ;
0075
0076 DataSet(N_Runs).data = data(energys == energy(j),:);
0077 % DataSet(N_Runs).data = DataSet(N_Runs).data .* (DataSet(N_Runs).data>zero_sup);
0078
0079 total_E = sum( DataSet(N_Runs).data* InitConst', 2) ;
0080
0081 DataSet(N_Runs).data = DataSet(N_Runs).data(total_E > 1, :);
0082 DataSet(N_Runs).accept = ones(size(data, 1), 1);
0083
0084 end
0085
0086 end
0087
0088 %%
0089 DrawDataSet(DataSet, InitConst, 'Inputs');
0090
0091 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0092
0093
0094
0095
0096 %%
0097
0098 figure('name',['DrawDataSet_EnergyFraction'],'PaperPositionMode','auto', ...
0099 'position',[100,0,1800,900]) ;
0100
0101 SumFraction = [];
0102
0103
0104 for i = 1:N_Runs
0105 subplot(1,2,i);
0106 total_E = ones(size(DataSet(i).data)) .*DataSet(i).E;
0107 Fraction = DataSet(i).data ./ total_E;
0108 MeanFraction = mean(Fraction, 1);
0109 MeanFraction = reshape(MeanFraction, 8, 8);
0110
0111 SumFraction = [SumFraction; Fraction];
0112
0113 imagesc(0:7, 0:7, MeanFraction);
0114 colorbar
0115 set(gca,'YDir','normal')
0116
0117 title(sprintf('%s, E = %.1f GeV', DataSet(i).FileID, DataSet(i).E));
0118 xlabel('Column ID');
0119 ylabel('Row ID');
0120
0121 end
0122
0123 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0124 %%
0125
0126
0127
0128 figure('name',['DrawDataSet_EnergyDistribution'],'PaperPositionMode','auto', ...
0129 'position',[100,0,1900,1400]) ;
0130
0131 for col = 1:8
0132 for row = 1:8
0133 idx = col + (row - 1)*8;
0134
0135 subplot(8,8,idx);
0136
0137 histogram( DataSet(i).data(:,idx),[0.:.1:16]);
0138
0139 set(gca,'XLim',[0,16])
0140 set(gca,'YScale','log')
0141 title(sprintf('Row%d Col%d', row-1, col-1));
0142
0143 end
0144 end
0145
0146 %%
0147
0148 figure('name',['DrawDataSet_EnergyFractionSum'],'PaperPositionMode','auto', ...
0149 'position',[100,0,1900,400]) ;
0150
0151 subplot(1,3,1);
0152
0153 MeanFraction = mean(SumFraction, 1);
0154 MeanFraction = reshape(MeanFraction, 8, 8);
0155
0156 imagesc(0:7, 0:7, MeanFraction);
0157 colorbar
0158 set(gca,'YDir','normal');
0159
0160 title(sprintf('Sum all data'));
0161 xlabel('Column ID');
0162 ylabel('Row ID');
0163
0164 subplot(1,3,2);
0165
0166 hist(reshape(MeanFraction,size(MeanFraction,1) * size(MeanFraction,2),1),100);
0167 title('fraction of shower energy carried by each tower');
0168 xlabel('MeanFraction');
0169
0170 subplot(1,3,3);
0171
0172 low_tower_IDs =reshape( MeanFraction<0.001, 1, Ndata);
0173 % low_tower_IDs =reshape( MeanFraction<0.0001, 1, Ndata);
0174
0175 imagesc(0:7, 0:7, reshape(low_tower_IDs, 8, 8));
0176 colorbar
0177 set(gca,'YDir','normal');
0178
0179 title(sprintf('Low hit towers'));
0180 xlabel('Column ID');
0181 ylabel('Row ID');
0182
0183 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0184
0185 % return
0186 %%
0187
0188 % InitConst_RunScale = [InitConst ones(1, N_Runs)];
0189 InitConst_RunScale = [InitConst ];
0190 % InitConst_RunScale = [InitConst ];
0191
0192 data_selection(InitConst_RunScale,10);
0193 disp(object_function(InitConst_RunScale));
0194
0195 % disp(object_function(InitConst_RunScale, 2));
0196 % disp(object_function(InitConst_RunScale, 1));
0197
0198 %%
0199
0200 x = InitConst_RunScale;
0201
0202 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0203
0204 data_selection(x,8);
0205 disp(object_function(x));
0206
0207 x = fminsearch(@(x) object_function(x), x,...
0208 options);
0209
0210 %%
0211 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0212 %
0213 data_selection(x,4);
0214 disp(object_function(x));
0215
0216 x = fminsearch(@(x) object_function(x), x,...
0217 options);
0218
0219
0220 %%
0221 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0222 %
0223 data_selection(x,2);
0224 disp(object_function(x));
0225
0226 x = fminsearch(@(x) object_function(x), x,...
0227 options);
0228
0229
0230 %%
0231
0232 data_selection(x,2);
0233
0234 calib_const = x(1:Ndata);
0235 % E_scale = x((Ndata+1):(Ndata + N_Runs));
0236 E_scale = ones(N_Runs);
0237
0238
0239 figure('name',['CalibConst'],'PaperPositionMode','auto', ...
0240 'position',[100,0,1800,400]) ;
0241
0242 subplot(1,3,1);
0243
0244 plot(calib_const);
0245 title(sprintf('Calibration constant'));
0246 xlabel('Col * 8 + Row');
0247 ylabel('Calibration New / Old');
0248
0249 subplot(1,3,2);
0250
0251 imagesc(0:7, 0:7, reshape(calib_const, 8, 8));
0252 colorbar
0253 set(gca,'YDir','normal')
0254
0255 title(sprintf('Calibration constant, New / Old'));
0256 xlabel('Column ID');
0257 ylabel('Row ID');
0258
0259
0260 subplot(1,3,3);
0261
0262 plot(E_scale,'x');
0263 title(sprintf('Energy scale constant, mean = %.2f', mean(E_scale)));
0264 xlabel('Run ID');
0265 ylabel('Energy scale New / Old');
0266 set(gca,'YLim',[0.5,1.5])
0267 grid on
0268
0269
0270 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0271 %%
0272
0273 % figure('name',['CalibConstVSModuleDensity'],'PaperPositionMode','auto', ...
0274 % 'position',[100,0,1300,1000]) ;
0275 %
0276 % ModuleDensity =[
0277 % 10.19 8.47 8.96 9.77
0278 % 9.74 9.96 10.00 9.87
0279 % 9.25 9.83 9.32 10.10
0280 % 9.59 9.39 10.06 9.62
0281 % 9.48 9.5 9.34 9.33
0282 % 9.43 9.34 9.39 9.55
0283 % 9.21 9.55 9.3 9.3
0284 % 9.5 9.6 9.3 9.24
0285 % ];
0286 %
0287 % ModuleDensity = ModuleDensity(8:-1:1,:);
0288 %
0289 % ModuleDensity = [ModuleDensity(:,1) ModuleDensity(:,1) ModuleDensity(:,2) ModuleDensity(:,2) ModuleDensity(:,3) ModuleDensity(:,3) ModuleDensity(:,4) ModuleDensity(:,4)];
0290 %
0291 % subplot(2,2,1);
0292 % imagesc(0:7, 0:7,ModuleDensity);
0293 % colorbar
0294 % set(gca,'YDir','normal')
0295 %
0296 % title(sprintf('Module Density (g/cm^3)'));
0297 % xlabel('Column ID');
0298 % ylabel('Row ID');
0299 %
0300 %
0301 % subplot(2,2,2);
0302 %
0303 % imagesc(0:7, 0:7, reshape(calib_const, 8, 8));
0304 % colorbar
0305 % set(gca,'YDir','normal');
0306 %
0307 % title(sprintf('Calibration constant, New / Old'));
0308 % xlabel('Column ID');
0309 % ylabel('Row ID');
0310 %
0311 %
0312 % subplot(2,2,3);
0313 %
0314 % [calib_const_col, calib_const_row]= meshgrid(0:7,0:7);
0315 %
0316 % calib_const_col = reshape(calib_const_col,1,64);
0317 % calib_const_row = reshape(calib_const_row,1,64);
0318 % % IDs = ~low_tower_IDs;
0319 % IDs = (calib_const_col >=3 & calib_const_col<=5 & calib_const_row>=4 & calib_const_row<=6);
0320 %
0321 % dens = reshape(ModuleDensity,1, 8* 8);
0322 % dens = dens(IDs);
0323 %
0324 % plot(calib_const(IDs),dens, 'o');
0325 % title(sprintf('THP 3x3, Calibration adjustment VS density'));
0326 % xlabel('Calibration constant, New / Old');
0327 % ylabel('Module Density');
0328 %
0329 % subplot(2,2,4);
0330 %
0331 %
0332 % IDs = ~low_tower_IDs;
0333 % % IDs = (calib_const_col >=3 & calib_const_col<=5 & calib_const_row>=4 & calib_const_row<=6);
0334 %
0335 % dens = reshape(ModuleDensity,1, 8* 8);
0336 % dens = dens(IDs);
0337 %
0338 % plot(calib_const(IDs),dens, 'o');
0339 % title(sprintf('All calibrated modules, Calibration adjustment VS density'));
0340 % xlabel('Calibration constant, New / Old');
0341 % ylabel('Module Density');
0342 %
0343 % SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0344
0345
0346
0347 %%
0348
0349 DrawDataSet(DataSet,calib_const,'Optimized');
0350 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0351
0352 %% Save out
0353 %
0354 % for i = 1:N_Runs
0355 % filename = sprintf('%s/calirbated_%d_%.0f.dat',DataFolder, DataSet(i).FileID, DataSet(i).E );
0356 %
0357 %
0358 % total_E = sum( DataSet(i).data*InitConst', 2) ;
0359 % calib_total_E = sum( DataSet(i).data* calib_const', 2) ;
0360 %
0361 % dlmwrite(filename,[total_E calib_total_E]);
0362 % end
0363
0364 save([DataFolder 'fit.mat']);
0365 % save('goodfit.mat');
0366
0367 %%
0368
0369 [calib_const_col, calib_const_row]= meshgrid(0:7,0:7);
0370
0371 A = [reshape(calib_const_col,1,64); reshape(calib_const_row,1,64); calib_const];
0372
0373
0374 fileID = fopen([DataFolder 'ShowerCalibFit_CablibConst.dat'],'w');
0375 % fprintf(fileID,'%d - %d \n',);
0376
0377 fprintf(fileID,'%d\t%d\t%f\n',A);
0378
0379
0380 fclose(fileID);
0381
0382 disp([DataFolder 'ShowerCalibFit_CablibConst.dat'])
0383
0384 %%
0385
0386 figure('name',['DrawDataSet_EnergyDistribution_AfterCalibration'],'PaperPositionMode','auto', ...
0387 'position',[100,0,1900,1400]) ;
0388
0389 for col = 1:8
0390 for row = 1:8
0391 idx = col + (row - 1)*8;
0392
0393 subplot(8,8,idx);
0394
0395 histogram( DataSet(i).data(:,idx) * calib_const(idx),[0.:.1:16]);
0396
0397 set(gca,'XLim',[0,16]);
0398 set(gca,'YScale','log');
0399 title(sprintf('Row%d Col%d', row-1, col-1));
0400
0401 end
0402 end
0403
