Warning, /analysis/Prototype3/EMCal/ShowerCalib/Fit/Proto3ShowerCalibFit.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 = {'3rd_positionscan'};
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(1,1,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.01, 1, Ndata);
0136 % low_tower_IDs =reshape( MeanFraction<0.0001, 1, Ndata);
0137
0138 imagesc(0:7, 0:7, reshape(low_tower_IDs, 8, 8));
0139 colorbar
0140 set(gca,'YDir','normal');
0141
0142 title(sprintf('Low hit towers'));
0143 xlabel('Column ID');
0144 ylabel('Row ID');
0145
0146 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0147
0148 % return
0149 %%
0150
0151 % InitConst_RunScale = [InitConst ones(1, N_Runs)];
0152 InitConst_RunScale = [InitConst ];
0153 % InitConst_RunScale = [InitConst ];
0154
0155 data_selection(InitConst_RunScale,10);
0156 disp(object_function(InitConst_RunScale));
0157
0158 % disp(object_function(InitConst_RunScale, 2));
0159 % disp(object_function(InitConst_RunScale, 1));
0160
0161 %%
0162
0163 x = InitConst_RunScale;
0164
0165 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0166
0167 data_selection(x,8);
0168 disp(object_function(x));
0169
0170 x = fminsearch(@(x) object_function(x), x,...
0171 options);
0172
0173 %%
0174 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0175 %
0176 data_selection(x,4);
0177 disp(object_function(x));
0178
0179 x = fminsearch(@(x) object_function(x), x,...
0180 options);
0181
0182
0183 %%
0184 options = optimset('Display','iter','TolX',1, 'MaxFunEvals', 100000,'MaxIter',40000,'PlotFcns',@optimplotfval );
0185 %
0186 data_selection(x,2);
0187 disp(object_function(x));
0188
0189 x = fminsearch(@(x) object_function(x), x,...
0190 options);
0191
0192
0193 %%
0194
0195 data_selection(x,2);
0196
0197 calib_const = x(1:Ndata);
0198 % E_scale = x((Ndata+1):(Ndata + N_Runs));
0199 E_scale = ones(N_Runs);
0200
0201
0202 figure('name',['CalibConst'],'PaperPositionMode','auto', ...
0203 'position',[100,0,1800,400]) ;
0204
0205 subplot(1,3,1);
0206
0207 plot(calib_const);
0208 title(sprintf('Calibration constant'));
0209 xlabel('Col * 8 + Row');
0210 ylabel('Calibration New / Old');
0211
0212 subplot(1,3,2);
0213
0214 imagesc(0:7, 0:7, reshape(calib_const, 8, 8));
0215 colorbar
0216 set(gca,'YDir','normal')
0217
0218 title(sprintf('Calibration constant, New / Old'));
0219 xlabel('Column ID');
0220 ylabel('Row ID');
0221
0222
0223 subplot(1,3,3);
0224
0225 plot(E_scale,'x');
0226 title(sprintf('Energy scale constant, mean = %.2f', mean(E_scale)));
0227 xlabel('Run ID');
0228 ylabel('Energy scale New / Old');
0229 set(gca,'YLim',[0.5,1.5])
0230 grid on
0231
0232
0233 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0234 %%
0235
0236 % figure('name',['CalibConstVSModuleDensity'],'PaperPositionMode','auto', ...
0237 % 'position',[100,0,1300,1000]) ;
0238 %
0239 % ModuleDensity =[
0240 % 10.19 8.47 8.96 9.77
0241 % 9.74 9.96 10.00 9.87
0242 % 9.25 9.83 9.32 10.10
0243 % 9.59 9.39 10.06 9.62
0244 % 9.48 9.5 9.34 9.33
0245 % 9.43 9.34 9.39 9.55
0246 % 9.21 9.55 9.3 9.3
0247 % 9.5 9.6 9.3 9.24
0248 % ];
0249 %
0250 % ModuleDensity = ModuleDensity(8:-1:1,:);
0251 %
0252 % ModuleDensity = [ModuleDensity(:,1) ModuleDensity(:,1) ModuleDensity(:,2) ModuleDensity(:,2) ModuleDensity(:,3) ModuleDensity(:,3) ModuleDensity(:,4) ModuleDensity(:,4)];
0253 %
0254 % subplot(2,2,1);
0255 % imagesc(0:7, 0:7,ModuleDensity);
0256 % colorbar
0257 % set(gca,'YDir','normal')
0258 %
0259 % title(sprintf('Module Density (g/cm^3)'));
0260 % xlabel('Column ID');
0261 % ylabel('Row ID');
0262 %
0263 %
0264 % subplot(2,2,2);
0265 %
0266 % imagesc(0:7, 0:7, reshape(calib_const, 8, 8));
0267 % colorbar
0268 % set(gca,'YDir','normal');
0269 %
0270 % title(sprintf('Calibration constant, New / Old'));
0271 % xlabel('Column ID');
0272 % ylabel('Row ID');
0273 %
0274 %
0275 % subplot(2,2,3);
0276 %
0277 % [calib_const_col, calib_const_row]= meshgrid(0:7,0:7);
0278 %
0279 % calib_const_col = reshape(calib_const_col,1,64);
0280 % calib_const_row = reshape(calib_const_row,1,64);
0281 % % IDs = ~low_tower_IDs;
0282 % IDs = (calib_const_col >=3 & calib_const_col<=5 & calib_const_row>=4 & calib_const_row<=6);
0283 %
0284 % dens = reshape(ModuleDensity,1, 8* 8);
0285 % dens = dens(IDs);
0286 %
0287 % plot(calib_const(IDs),dens, 'o');
0288 % title(sprintf('THP 3x3, Calibration adjustment VS density'));
0289 % xlabel('Calibration constant, New / Old');
0290 % ylabel('Module Density');
0291 %
0292 % subplot(2,2,4);
0293 %
0294 %
0295 % IDs = ~low_tower_IDs;
0296 % % IDs = (calib_const_col >=3 & calib_const_col<=5 & calib_const_row>=4 & calib_const_row<=6);
0297 %
0298 % dens = reshape(ModuleDensity,1, 8* 8);
0299 % dens = dens(IDs);
0300 %
0301 % plot(calib_const(IDs),dens, 'o');
0302 % title(sprintf('All calibrated modules, Calibration adjustment VS density'));
0303 % xlabel('Calibration constant, New / Old');
0304 % ylabel('Module Density');
0305 %
0306 % SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0307
0308
0309
0310 %%
0311
0312 DrawDataSet(DataSet,calib_const,'Optimized');
0313 SaveCanvas([DataFolder 'EnergyCalibFIt'],gcf);
0314
0315 %% Save out
0316 %
0317 % for i = 1:N_Runs
0318 % filename = sprintf('%s/calirbated_%d_%.0f.dat',DataFolder, DataSet(i).FileID, DataSet(i).E );
0319 %
0320 %
0321 % total_E = sum( DataSet(i).data*InitConst', 2) ;
0322 % calib_total_E = sum( DataSet(i).data* calib_const', 2) ;
0323 %
0324 % dlmwrite(filename,[total_E calib_total_E]);
0325 % end
0326
0327 save([DataFolder 'fit.mat']);
0328 % save('goodfit.mat');
0329
0330 %%
0331
0332 [calib_const_col, calib_const_row]= meshgrid(0:7,0:7);
0333
0334 A = [reshape(calib_const_col,1,64); reshape(calib_const_row,1,64); calib_const];
0335
0336 fileID = fopen([DataFolder 'ShowerCalibFit_CablibConst.dat'],'w');
0337 % fprintf(fileID,'%d - %d \n',);
0338
0339 fprintf(fileID,'%d\t%d\t%f\n',A);
0340
0341
0342 fclose(fileID);
0343
0344
0345
0346
