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