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); |
---|