[17] | 1 | function [QMS, WarningString] = quadplot(Input1, FigureHandle, sigmaBPM) |
---|
| 2 | %QUADPLOT - Plots quadrupole centering data |
---|
| 3 | % [QMS, WarningString] = quadplot(Input1, Handle, sigmaBPM) |
---|
| 4 | % |
---|
| 5 | % INPUTS |
---|
| 6 | % 1. Input1 can be a filename for the data or a QMS structure (see help quadcenter for details). |
---|
| 7 | % If empty or zero inputs, then a dialog box will be provided to select a file. |
---|
| 8 | % 2. Handle can be a figure handle or a vector of 4 axes handles |
---|
| 9 | % If Handle=0, no results are plotted |
---|
| 10 | % 3. Standard deviation of the BPMs (scalar if all BPMs are the same) |
---|
| 11 | % These should be in the data file, but this provides an override if not |
---|
| 12 | % found then the default is inf (ie, not used). |
---|
| 13 | % |
---|
| 14 | % OUTPUTS |
---|
| 15 | % 1. For details of the QMS data structure see help quadcenter |
---|
| 16 | % This function added: |
---|
| 17 | % QMS.Offset - Offset computed at each BPM |
---|
| 18 | % QMS.FitParameters - Fit parameter at each BPM |
---|
| 19 | % QMS.FitParametersStd - Sigma of the fit parameter at each BPM |
---|
| 20 | % QMS.BPMStd - BPM sigma at each BPM |
---|
| 21 | % QMS.OffsetSTDMontiCarlo - Monti Carlo estimate of the sigma of the offset (optional) |
---|
| 22 | % |
---|
| 23 | % 2. WarningString = string with warning message if you occurred |
---|
| 24 | % |
---|
| 25 | % See also quadcenter, quaderrors |
---|
| 26 | |
---|
| 27 | % |
---|
| 28 | %% Written by Greg Portmann |
---|
| 29 | |
---|
| 30 | % Laurent Out |
---|
| 31 | % |
---|
| 32 | % To Do: |
---|
| 33 | % 1. It wouldn't be to difficult to added a LS weight based on slope or even the ideal weighting of std(center). |
---|
| 34 | % I haven't done it yet because the BPM errors are usually roughly equal at most accelerators. |
---|
| 35 | |
---|
| 36 | |
---|
| 37 | % Remove BPM gain and roll before finding the center |
---|
| 38 | ModelCoordinates = 0; |
---|
| 39 | |
---|
| 40 | % Figure setup |
---|
| 41 | Buffer = .03; |
---|
| 42 | HeightBuffer = .08; |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | % Remove BPM if it's slope less than MinSlopeFraction * (the maximum slope) |
---|
| 46 | MinSlopeFraction = .25; |
---|
| 47 | |
---|
| 48 | % # of STD of the center calculation allowed |
---|
| 49 | CenterOutlierFactor = 1; |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | QMS = []; |
---|
| 53 | WarningString = ''; |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | % Inputs |
---|
| 57 | try |
---|
| 58 | if nargin == 0 |
---|
| 59 | [FileName, PathName] = uigetfile('*.mat', 'Select a Quadrupole Center File', [getfamilydata('Directory','DataRoot'), 'QMS', filesep]); |
---|
| 60 | if ~isstr(FileName) |
---|
| 61 | return |
---|
| 62 | else |
---|
| 63 | load([PathName,FileName]); |
---|
| 64 | end |
---|
| 65 | else |
---|
| 66 | if isempty(Input1) |
---|
| 67 | [FileName, PathName] = uigetfile('*.mat', 'Select a Quadrupole Center File'); |
---|
| 68 | %[FileName, PathName] = uigetfile('*.mat', 'Select a Quadrupole Center File', [getfamilydata('Directory','DataRoot'), 'QMS', filesep]); |
---|
| 69 | if ~isstr(FileName) |
---|
| 70 | return |
---|
| 71 | else |
---|
| 72 | load([PathName,FileName]); |
---|
| 73 | end |
---|
| 74 | elseif isstr(Input1) |
---|
| 75 | FileName = Input1; |
---|
| 76 | load(FileName); |
---|
| 77 | else |
---|
| 78 | QMS = Input1; |
---|
| 79 | FileName = []; |
---|
| 80 | end |
---|
| 81 | end |
---|
| 82 | catch |
---|
| 83 | error('Problem getting input data'); |
---|
| 84 | end |
---|
| 85 | if nargin < 2 |
---|
| 86 | FigureHandle = []; |
---|
| 87 | end |
---|
| 88 | |
---|
| 89 | [BPMelem1, iNotFound] = findrowindex(QMS.BPMDev, QMS.BPMDevList); |
---|
| 90 | if ~isempty(iNotFound) |
---|
| 91 | error('BPM at the quadrupole not found in the BPM device list'); |
---|
| 92 | end |
---|
| 93 | |
---|
| 94 | QuadraticFit = QMS.QuadraticFit; % 0 = linear fit, else quadratic fit |
---|
| 95 | |
---|
| 96 | % TO SEE |
---|
| 97 | % At SOLEIL noise is very small, constraint too large |
---|
| 98 | %OutlierFactor = QMS.OutlierFactor; % if abs(data - fit) > OutlierFactor*BPMsdt, then remove that BPM |
---|
| 99 | OutlierFactor = 1e-3; % if abs(data - fit) > OutlierFactor, then remove that BPM |
---|
| 100 | |
---|
| 101 | % Get BPM standard deviation |
---|
| 102 | if nargin < 3 |
---|
| 103 | % Get from the data file |
---|
| 104 | if isfield(QMS, 'BPMSTD') |
---|
| 105 | sigmaBPM = QMS.BPMSTD; |
---|
| 106 | else |
---|
| 107 | sigmaBPM = inf; |
---|
| 108 | end |
---|
| 109 | end |
---|
| 110 | if isempty(sigmaBPM) |
---|
| 111 | sigmaBPM = inf; |
---|
| 112 | end |
---|
| 113 | if isnan(sigmaBPM) | isinf(sigmaBPM) |
---|
| 114 | sigmaBPM = inf; |
---|
| 115 | fprintf(' WARNING: BPM standard deviation is unknown, hence there is no BPM outlier condition.\n'); |
---|
| 116 | end |
---|
| 117 | sigmaBPM = sigmaBPM(:); |
---|
| 118 | QMS.BPMSTD = sigmaBPM; |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | % Get figure handle |
---|
| 122 | if all(FigureHandle ~= 0) |
---|
| 123 | if isempty(FigureHandle) |
---|
| 124 | FigureHandle = figure; |
---|
| 125 | clf reset |
---|
| 126 | AxesHandles(1) = subplot(4,1,1); |
---|
| 127 | AxesHandles(2) = subplot(4,1,2); |
---|
| 128 | AxesHandles(3) = subplot(4,1,3); |
---|
| 129 | AxesHandles(4) = subplot(4,1,4); |
---|
| 130 | else |
---|
| 131 | if length(FigureHandle) == 1 |
---|
| 132 | FigureHandle = figure(FigureHandle); |
---|
| 133 | clf reset |
---|
| 134 | AxesHandles(1) = subplot(4,1,1); |
---|
| 135 | AxesHandles(2) = subplot(4,1,2); |
---|
| 136 | AxesHandles(3) = subplot(4,1,3); |
---|
| 137 | AxesHandles(4) = subplot(4,1,4); |
---|
| 138 | elseif length(FigureHandle) == 4 |
---|
| 139 | FigureHandle = figure; |
---|
| 140 | clf reset |
---|
| 141 | AxesHandles = FigureHandle; |
---|
| 142 | else |
---|
| 143 | error('Improper size of input FigureHandle'); |
---|
| 144 | end |
---|
| 145 | end |
---|
| 146 | end |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | % Change the coordinates to model X & Y |
---|
| 151 | if ModelCoordinates |
---|
| 152 | fprintf(' Changing to model coordinates (the final center should be "rotated" back to the raw BPM coordinates).\n'); |
---|
| 153 | Gx = getgain('BPMx', QMS.BPMDevList); % Unfortunately I don't have both families in the QMS structure??? |
---|
| 154 | Gy = getgain('BPMy', QMS.BPMDevList); |
---|
| 155 | Crunch = getcrunch(QMS.BPMFamily, QMS.BPMDevList); |
---|
| 156 | Roll = getroll(QMS.BPMFamily, QMS.BPMDevList); |
---|
| 157 | |
---|
| 158 | for i = 1:length(Gx) |
---|
| 159 | M = gcr2loco(Gx(i), Gy(i), Crunch(i), Roll(i)); |
---|
| 160 | invM = inv(M); |
---|
| 161 | |
---|
| 162 | tmp = invM * [QMS.x0(i,:); QMS.y0(i,:)]; |
---|
| 163 | QMS.x0(i,:) = tmp(1,:); |
---|
| 164 | QMS.y0(i,:) = tmp(2,:); |
---|
| 165 | |
---|
| 166 | tmp = invM * [QMS.x1(i,:); QMS.y1(i,:)]; |
---|
| 167 | QMS.x1(i,:) = tmp(1,:); |
---|
| 168 | QMS.y1(i,:) = tmp(2,:); |
---|
| 169 | |
---|
| 170 | tmp = invM * [QMS.x2(i,:); QMS.y2(i,:)]; |
---|
| 171 | QMS.x2(i,:) = tmp(1,:); |
---|
| 172 | QMS.y2(i,:) = tmp(2,:); |
---|
| 173 | |
---|
| 174 | tmp = invM * [QMS.xstart(i,:); QMS.ystart(i,:)]; |
---|
| 175 | QMS.xstart(i,:) = tmp(1,:); |
---|
| 176 | QMS.ystart(i,:) = tmp(2,:); |
---|
| 177 | end |
---|
| 178 | end |
---|
| 179 | |
---|
| 180 | |
---|
| 181 | % Select the x or y plane |
---|
| 182 | if QMS.QuadPlane == 1 |
---|
| 183 | x0 = QMS.x0; |
---|
| 184 | x1 = QMS.x1; |
---|
| 185 | x2 = QMS.x2; |
---|
| 186 | |
---|
| 187 | % Plot setup |
---|
| 188 | if all(FigureHandle ~= 0) |
---|
| 189 | set(FigureHandle,'units','normal','position',[Buffer .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]); |
---|
| 190 | end |
---|
| 191 | else |
---|
| 192 | x0 = QMS.y0; |
---|
| 193 | x1 = QMS.y1; |
---|
| 194 | x2 = QMS.y2; |
---|
| 195 | |
---|
| 196 | % Plot setup |
---|
| 197 | if all(FigureHandle ~= 0) |
---|
| 198 | set(FigureHandle,'units','normal','position',[.503 .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]); |
---|
| 199 | end |
---|
| 200 | end |
---|
| 201 | |
---|
| 202 | |
---|
| 203 | |
---|
| 204 | % Expand sigmaBPM is necessary |
---|
| 205 | if length(sigmaBPM) == 1 |
---|
| 206 | sigmaBPM = ones(size(x1,1),1) * sigmaBPM; |
---|
| 207 | end |
---|
| 208 | |
---|
| 209 | N = size(x1,2); |
---|
| 210 | |
---|
| 211 | |
---|
| 212 | % Fit verses the position at the BPM next to the quadrupole |
---|
| 213 | if strcmpi(QMS.ModulationMethod, 'sweep') |
---|
| 214 | % One dimensional sweep of the quadrupole |
---|
| 215 | %x = x1(BPMelem1,:)'; |
---|
| 216 | x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2; |
---|
| 217 | elseif strcmpi(QMS.ModulationMethod, 'bipolar') |
---|
| 218 | % Modulation of the quadrupole |
---|
| 219 | x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2; |
---|
| 220 | elseif strcmpi(QMS.ModulationMethod, 'unipolar') |
---|
| 221 | % Unipolar modulation of the quadrupole |
---|
| 222 | %x = x1(BPMelem1,:)'; |
---|
| 223 | x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2; |
---|
| 224 | end |
---|
| 225 | |
---|
| 226 | if isfield(QMS, 'Orbit0') |
---|
| 227 | BPMUnitsString = QMS.Orbit0.UnitsString; |
---|
| 228 | else |
---|
| 229 | BPMUnitsString = 'mm'; |
---|
| 230 | end |
---|
| 231 | |
---|
| 232 | % Figure #1 |
---|
| 233 | if all(FigureHandle ~= 0) |
---|
| 234 | |
---|
| 235 | % Plot horizontal raw data |
---|
| 236 | axes(AxesHandles(1)); |
---|
| 237 | plot((QMS.x2(BPMelem1,:)+QMS.x1(BPMelem1,:))/2, (QMS.x2-QMS.x1), '-x'); |
---|
| 238 | hold on; |
---|
| 239 | % just data with corrector value at Inominal - deltaI |
---|
| 240 | plot((QMS.x2(BPMelem1,1)+QMS.x1(BPMelem1,1))/2, (QMS.x2(:,1)-QMS.x1(:,1)), 'xb'); |
---|
| 241 | hold off; |
---|
| 242 | xlabel(sprintf('%s(%d,%d), raw values [%s]', 'BPMx', QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString)); |
---|
| 243 | ylabel({sprintf('\\Delta %s [%s]', 'BPMx', BPMUnitsString),'(raw data)'}); |
---|
| 244 | if isempty(FileName) |
---|
| 245 | title(sprintf('Center for %s(%d,%d) %s(%d,%d)', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev), 'interpreter', 'none'); |
---|
| 246 | else |
---|
| 247 | title(sprintf('Center for %s(%d,%d) %s(%d,%d) (%s)', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev, FileName), 'interpreter', 'none'); |
---|
| 248 | end |
---|
| 249 | grid on |
---|
| 250 | axis tight |
---|
| 251 | |
---|
| 252 | if QMS.QuadPlane == 1 |
---|
| 253 | Axes2Axis = axis; |
---|
| 254 | end |
---|
| 255 | |
---|
| 256 | % Plot vertical raw data |
---|
| 257 | axes(AxesHandles(2)); |
---|
| 258 | plot((QMS.y2(BPMelem1,:)+QMS.y1(BPMelem1,:))/2, (QMS.y2-QMS.y1), '-x'); |
---|
| 259 | hold on; |
---|
| 260 | plot((QMS.y2(BPMelem1,1)+QMS.y1(BPMelem1,1))/2, (QMS.y2(:,1)-QMS.y1(:,1)), 'xb'); |
---|
| 261 | hold off; |
---|
| 262 | %plot(x, x2-x1, '-x'); |
---|
| 263 | %plot(linspace(-DelHCM,DelHCM,3), x2-x1); |
---|
| 264 | xlabel(sprintf('%s(%d,%d), raw values [%s]', 'BPMy', QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString)); |
---|
| 265 | ylabel({sprintf('\\Delta %s [%s]', 'BPMy', BPMUnitsString),'(raw data)'}); |
---|
| 266 | grid on |
---|
| 267 | axis tight |
---|
| 268 | |
---|
| 269 | if QMS.QuadPlane == 2 |
---|
| 270 | Axes2Axis = axis; |
---|
| 271 | end |
---|
| 272 | end |
---|
| 273 | |
---|
| 274 | |
---|
| 275 | if isempty(FileName) |
---|
| 276 | fprintf(' Calculating the center of %s(%d,%d) using %s(%d,%d)\n', QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev); |
---|
| 277 | else |
---|
| 278 | fprintf(' Calculating the center of %s(%d,%d) using %s(%d,%d) (Data file: %s)\n', QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev, FileName); |
---|
| 279 | end |
---|
| 280 | fprintf(' Quadrupole modulation delta = %.3f amps, %s(%d,%d) max step = %.3f amps (%s)\n', QMS.QuadDelta, QMS.CorrFamily, QMS.CorrDevList(1,:), QMS.CorrDelta, QMS.ModulationMethod); |
---|
| 281 | |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | % Least squares fit |
---|
| 285 | merit = x2 - x1; |
---|
| 286 | if QuadraticFit |
---|
| 287 | X = [ones(size(x)) x x.^2]; |
---|
| 288 | fprintf(' %d point parabolic least squares fit\n', N); |
---|
| 289 | else |
---|
| 290 | X = [ones(size(x)) x]; |
---|
| 291 | fprintf(' %d point linear least squares fit\n', N); |
---|
| 292 | end |
---|
| 293 | |
---|
| 294 | |
---|
| 295 | % Axes #3 |
---|
| 296 | if all(FigureHandle ~= 0) |
---|
| 297 | axes(AxesHandles(3)); |
---|
| 298 | xx = linspace(x(1), x(end), 200); |
---|
| 299 | end |
---|
| 300 | |
---|
| 301 | iBPMOutlier = []; |
---|
| 302 | invXX = inv(X'*X); |
---|
| 303 | invXX_X = invXX * X'; |
---|
| 304 | |
---|
| 305 | for i = 1:size(x1,1) |
---|
| 306 | % least-square fit: m = slope and b = Y-intercept |
---|
| 307 | b = invXX_X * merit(i,:)'; |
---|
| 308 | bhat(i,:) = b'; |
---|
| 309 | |
---|
| 310 | % Should equal |
---|
| 311 | %b = X \merit(i,:)'; |
---|
| 312 | %bhat1(i,:) = b'; |
---|
| 313 | |
---|
| 314 | % Standard deviation |
---|
| 315 | bstd = sigmaBPM(i) * invXX; |
---|
| 316 | bhatstd(i,:) = diag(bstd)'; % hopefully cross-correlation terms are small |
---|
| 317 | |
---|
| 318 | if all(FigureHandle ~= 0) |
---|
| 319 | if QuadraticFit |
---|
| 320 | y = b(3)*xx.^2 + b(2)*xx + b(1); |
---|
| 321 | else |
---|
| 322 | y = b(2)*xx + b(1); |
---|
| 323 | end |
---|
| 324 | % plot(xx, y); hold on |
---|
| 325 | end |
---|
| 326 | |
---|
| 327 | % Outlier condition: remove if the error between the fit and the data is greater than OutlierFactor |
---|
| 328 | if QuadraticFit |
---|
| 329 | y = b(3)*x.^2 + b(2)*x + b(1); |
---|
| 330 | else |
---|
| 331 | y = b(2)*x + b(1); |
---|
| 332 | end |
---|
| 333 | if max(abs(y - merit(i,:)')) > OutlierFactor; % OutlierFactor was absolute max(abs(y - merit(i,:)')) > OutlierFactor |
---|
| 334 | %if max(abs(y - merit(i,:)')) > OutlierFactor * sigmaBPM(i) %Laurent |
---|
| 335 | iBPMOutlier = [iBPMOutlier;i]; |
---|
| 336 | end |
---|
| 337 | |
---|
| 338 | if QuadraticFit |
---|
| 339 | % Quadratic fit |
---|
| 340 | if b(2) > 0 |
---|
| 341 | offset(i,1) = (-b(2) + sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3)); |
---|
| 342 | else |
---|
| 343 | offset(i,1) = (-b(2) - sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3)); |
---|
| 344 | end |
---|
| 345 | if ~isreal(offset(i,1)) |
---|
| 346 | % (b^2-4ac) can be negative but it will only happen if the slope is very small. The offset |
---|
| 347 | % should just get thrown out later as an outlier but change the solution to the minimum of the parabola. |
---|
| 348 | offset(i,1) = -b(2) / b(1) / 2; |
---|
| 349 | end |
---|
| 350 | else |
---|
| 351 | % Linear fit |
---|
| 352 | offset(i,1) = -b(1)/b(2); |
---|
| 353 | end |
---|
| 354 | end |
---|
| 355 | |
---|
| 356 | QMS.Offset = offset; |
---|
| 357 | QMS.FitParameters = bhat; |
---|
| 358 | QMS.FitParametersStd = bhatstd; |
---|
| 359 | QMS.BPMStd = sigmaBPM; |
---|
| 360 | |
---|
| 361 | |
---|
| 362 | % % Label axes #2 |
---|
| 363 | % if all(FigureHandle ~= 0) |
---|
| 364 | % xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString)); |
---|
| 365 | % ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString)); |
---|
| 366 | % grid on |
---|
| 367 | % axis tight |
---|
| 368 | % end |
---|
| 369 | |
---|
| 370 | |
---|
| 371 | % Compute offset for different conditions |
---|
| 372 | fprintf(' %d total BPMs\n', length(offset)); |
---|
| 373 | fprintf(' BPMs are removed for 1. Bad Status, 2. BPM Outlier, 3. Small Slope, or 4. Center Outlier\n'); |
---|
| 374 | %m1 = mean(offset); |
---|
| 375 | %s1 = std(offset); |
---|
| 376 | %fprintf(' 0. Mean = %.5f %s, STD = %.5f %s, all %d BPMs\n', m1, BPMUnitsString, s1, BPMUnitsString, length(offset)); |
---|
| 377 | |
---|
| 378 | |
---|
| 379 | % Remove bad Status BPMs |
---|
| 380 | iStatus = find(QMS.BPMStatus==0); |
---|
| 381 | iBad = iStatus; |
---|
| 382 | if length(iBad) == length(offset) |
---|
| 383 | error('All the BPMs have bad status'); |
---|
| 384 | end |
---|
| 385 | offset1 = offset; |
---|
| 386 | offset1(iBad) = []; |
---|
| 387 | m2 = mean(offset1); |
---|
| 388 | s2 = std(offset1); |
---|
| 389 | fprintf(' 1. Mean = %+.5f %s, STD = %.5f %s, %2d points with bad status\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad)); |
---|
| 390 | |
---|
| 391 | % Remove bad Status + Outliers |
---|
| 392 | iBad = unique([iStatus; iBPMOutlier]); |
---|
| 393 | if length(iBad) == length(offset) |
---|
| 394 | error('All BPMs either have bad status or failed the BPM outlier condition'); |
---|
| 395 | end |
---|
| 396 | offset1 = offset; |
---|
| 397 | offset1(iBad) = []; |
---|
| 398 | m2 = mean(offset1); |
---|
| 399 | s2 = std(offset1); |
---|
| 400 | %fprintf(' 2. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1 and abs(fit - measured data) > %.2f std(BPM) (BPM outlier)\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), OutlierFactor); |
---|
| 401 | fprintf(' 2. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1 and abs(fit - measured data) > %.4f (BPM outlier)\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), OutlierFactor); |
---|
| 402 | |
---|
| 403 | |
---|
| 404 | % Remove bad Status + Small slopes |
---|
| 405 | %iSlope = find(abs(bhat(:,2)) < max(abs(bhat(:,2)))*MinSlopeFraction); |
---|
| 406 | |
---|
| 407 | % Look for slope outliers |
---|
| 408 | Slopes = abs(bhat(:,2)); |
---|
| 409 | [Slopes, i] = sort(Slopes); |
---|
| 410 | Slopes = Slopes(round(end/2):end); % remove the first half |
---|
| 411 | if length(Slopes) > 5 |
---|
| 412 | SlopesMax = Slopes(end-4); |
---|
| 413 | else |
---|
| 414 | SlopesMax = Slopes(end); |
---|
| 415 | end |
---|
| 416 | %i = find(abs(Slopes-mean(Slopes)) > 2 * std(Slopes)); |
---|
| 417 | %Slopes(i) = []; |
---|
| 418 | |
---|
| 419 | iSlope = find(abs(bhat(:,2)) < SlopesMax * MinSlopeFraction); |
---|
| 420 | iBad = unique([iStatus; iBPMOutlier; iSlope]); |
---|
| 421 | if length(iBad) == length(offset) |
---|
| 422 | error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition'); |
---|
| 423 | end |
---|
| 424 | offset1 = offset; |
---|
| 425 | offset1(iBad) = []; |
---|
| 426 | m2 = mean(offset1); |
---|
| 427 | s2 = std(offset1); |
---|
| 428 | fprintf(' 3. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1, 2, and slope < %.2f max(slope)\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), MinSlopeFraction); |
---|
| 429 | |
---|
| 430 | |
---|
| 431 | % Offset outlier offsets-mean(offsets) greater than 1 std |
---|
| 432 | itotal = (1:length(offset))'; |
---|
| 433 | iok = itotal; |
---|
| 434 | |
---|
| 435 | offset1 = offset; |
---|
| 436 | offset1(iBad) = []; |
---|
| 437 | iok(iBad) = []; |
---|
| 438 | |
---|
| 439 | i = find(abs(offset1-mean(offset1)) > CenterOutlierFactor * std(offset1)); |
---|
| 440 | iCenterOutlier = iok(i); |
---|
| 441 | iBad = unique([iBad; iCenterOutlier]); |
---|
| 442 | if length(iBad) == length(offset) |
---|
| 443 | error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition, , or failed the center outlier condition'); |
---|
| 444 | end |
---|
| 445 | offset1(i) = []; |
---|
| 446 | iok(i) = []; |
---|
| 447 | |
---|
| 448 | m2 = mean(offset1); |
---|
| 449 | s2 = std(offset1); |
---|
| 450 | QMS.GoodIndex = iok; |
---|
| 451 | fprintf(' 4. Mean = %+.5f %s, STD = %.5f %s, %2d points with condition 1, 2, 3, and abs(center-mean(center)) > %.2f std(center) (Center outlier)\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad), CenterOutlierFactor); |
---|
| 452 | |
---|
| 453 | |
---|
| 454 | NN = length(offset); |
---|
| 455 | |
---|
| 456 | % Axes #4 |
---|
| 457 | if all(FigureHandle ~= 0) |
---|
| 458 | axes(AxesHandles(4)); |
---|
| 459 | [xx1,yy1]=stairs(1:NN,offset); |
---|
| 460 | offset1 = offset; |
---|
| 461 | offset1(iBad) = NaN*ones(length(iBad),1); |
---|
| 462 | [xx2, yy2] = stairs(1:NN, offset1); |
---|
| 463 | plot(xx1,yy1,'r', xx2,yy2,'b'); |
---|
| 464 | xlabel('BPM Number'); |
---|
| 465 | ylabel(sprintf('%s Center [%s]', QMS.BPMFamily, BPMUnitsString)); |
---|
| 466 | grid on |
---|
| 467 | axis tight |
---|
| 468 | %xaxis([0 NN+1]); |
---|
| 469 | axis([0 NN+1 min(offset1)-.1e-3 max(offset1)+.1e-3]); |
---|
| 470 | end |
---|
| 471 | |
---|
| 472 | |
---|
| 473 | % Axes #3 |
---|
| 474 | |
---|
| 475 | if all(FigureHandle ~= 0) |
---|
| 476 | if 0 |
---|
| 477 | % Plot red line over the bad lines |
---|
| 478 | axes(AxesHandles(3)); |
---|
| 479 | for j = 1:length(iBad) |
---|
| 480 | i = iBad(j); |
---|
| 481 | if QuadraticFit |
---|
| 482 | y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1); |
---|
| 483 | else |
---|
| 484 | y = bhat(i,2)*xx + bhat(i,1); |
---|
| 485 | end |
---|
| 486 | plot(xx, y,'r'); |
---|
| 487 | end |
---|
| 488 | hold off |
---|
| 489 | axis tight |
---|
| 490 | else |
---|
| 491 | % Only plot the good data |
---|
| 492 | axes(AxesHandles(3)); |
---|
| 493 | yy = []; |
---|
| 494 | for i = 1:size(x1,1) |
---|
| 495 | if ~any(i == iBad) |
---|
| 496 | if QuadraticFit |
---|
| 497 | y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1); |
---|
| 498 | else |
---|
| 499 | y = bhat(i,2)*xx + bhat(i,1); |
---|
| 500 | end |
---|
| 501 | yy = [yy;y]; |
---|
| 502 | end |
---|
| 503 | end |
---|
| 504 | plot(xx, yy); |
---|
| 505 | hold off |
---|
| 506 | grid on |
---|
| 507 | axis tight |
---|
| 508 | end |
---|
| 509 | xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString)); |
---|
| 510 | ylabel({sprintf('\\Delta %s [%s]', QMS.BPMFamily, BPMUnitsString),'(LS fit)'}); |
---|
| 511 | grid on |
---|
| 512 | axis tight |
---|
| 513 | |
---|
| 514 | xaxis(Axes2Axis(1:2)); |
---|
| 515 | end |
---|
| 516 | |
---|
| 517 | |
---|
| 518 | |
---|
| 519 | if ~isempty(iStatus) |
---|
| 520 | if ~isempty(find(iStatus==BPMelem1)) |
---|
| 521 | fprintf(' WARNING: BPM(%d,%d) has a bad status\n', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 522 | WarningString = sprintf('BPM(%d,%d) has a bad status', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 523 | end |
---|
| 524 | end |
---|
| 525 | if ~isempty(iBPMOutlier) |
---|
| 526 | if ~isempty(find(iBPMOutlier==BPMelem1)) |
---|
| 527 | fprintf(' WARNING: BPM(%d,%d) removed due to outlier (based on std(BPM))\n', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 528 | WarningString = sprintf('BPM(%d,%d) removed due to outlier (based on std(BPM))', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 529 | end |
---|
| 530 | end |
---|
| 531 | if ~isempty(iSlope) |
---|
| 532 | if ~isempty(find(iSlope==BPMelem1)) |
---|
| 533 | fprintf(' WARNING: BPM(%d,%d) slope is too small\n', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 534 | WarningString = sprintf('BPM(%d,%d) slope is too small', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 535 | end |
---|
| 536 | end |
---|
| 537 | if ~isempty(iCenterOutlier) |
---|
| 538 | if ~isempty(find(iCenterOutlier==BPMelem1)) |
---|
| 539 | fprintf(' WARNING: BPM(%d,%d) removed due to outlier (based on all the centers)\n', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 540 | WarningString = sprintf('BPM(%d,%d) ', QMS.BPMDev(1), QMS.BPMDev(2)); |
---|
| 541 | end |
---|
| 542 | end |
---|
| 543 | |
---|
| 544 | |
---|
| 545 | % % Axes #3 |
---|
| 546 | % if all(FigureHandle ~= 0) |
---|
| 547 | % axes(AxesHandles(4)); |
---|
| 548 | % iii=1:NN; |
---|
| 549 | % iii(iBad)=[]; |
---|
| 550 | % for j = 1:length(iii) |
---|
| 551 | % i = iii(j); |
---|
| 552 | % |
---|
| 553 | % if 1 |
---|
| 554 | % % Plot fit |
---|
| 555 | % if QuadraticFit |
---|
| 556 | % y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1); |
---|
| 557 | % else |
---|
| 558 | % y = bhat(i,2)*xx + bhat(i,1); |
---|
| 559 | % end |
---|
| 560 | % if all(FigureHandle ~= 0) |
---|
| 561 | % plot(xx, y,'b'); hold on |
---|
| 562 | % end |
---|
| 563 | % else |
---|
| 564 | % % Plot error in fit |
---|
| 565 | % if QuadraticFit |
---|
| 566 | % y = bhat(i,3)*x.^2 + bhat(i,2)*x + bhat(i,1); |
---|
| 567 | % else |
---|
| 568 | % y = bhat(i,2)*x + bhat(i,1); |
---|
| 569 | % end |
---|
| 570 | % plot(x, y - merit(i,:)','b'); hold on |
---|
| 571 | % end |
---|
| 572 | % end |
---|
| 573 | % hold off |
---|
| 574 | % xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString)); |
---|
| 575 | % ylabel(sprintf('Final %s Merit Fun [%s]', QMS.BPMFamily, BPMUnitsString)); |
---|
| 576 | % grid on |
---|
| 577 | % axis tight |
---|
| 578 | % orient tall |
---|
| 579 | % end |
---|
| 580 | |
---|
| 581 | |
---|
| 582 | if ~isempty(QMS.OldCenter) |
---|
| 583 | fprintf(' Starting Offset %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, QMS.OldCenter, BPMUnitsString); |
---|
| 584 | end |
---|
| 585 | fprintf(' New Offset %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, m2, BPMUnitsString); |
---|
| 586 | |
---|
| 587 | QMS.Center = m2; |
---|
| 588 | QMS.CenterSTD = s2; |
---|
| 589 | |
---|
| 590 | |
---|
| 591 | if all(FigureHandle ~= 0) |
---|
| 592 | addlabel(1, 0, datestr(QMS.TimeStamp)); |
---|
| 593 | addlabel(0, 0, sprintf('Offset %s(%d,%d)=%f, {\\Delta}%s(%d,%d)=%f, {\\Delta}%s(%d,%d)=%f', QMS.BPMFamily, QMS.BPMDev, QMS.Center, QMS.QuadFamily, QMS.QuadDev, QMS.QuadDelta, QMS.CorrFamily, QMS.CorrDevList(1,:), QMS.CorrDelta(1)),10); |
---|
| 594 | end |
---|
| 595 | |
---|
| 596 | fprintf('\n'); |
---|
| 597 | |
---|
| 598 | |
---|
| 599 | % Get and plot errors |
---|
| 600 | if all(FigureHandle ~= 0) |
---|
| 601 | QMS = quaderrors(QMS, gcf+1); |
---|
| 602 | else |
---|
| 603 | % This is a cluge for now |
---|
| 604 | QMS = quaderrors(QMS, 0); |
---|
| 605 | end |
---|
| 606 | |
---|
| 607 | |
---|
| 608 | %QMS = orderfields(QMS); |
---|