quadplotphysics

QUADPLOTPHYSICS - Plots quadrupole centering data in physics units
  [QMS, WarningString] = quadplot(Input1, Handle, sigmaBPM)

  INPUTS
  1. Input1 can be a filename for the data or a QMS structure (see help quadcenter for details).
     If empty or zero inputs, then a dialog box will be provided to select a file.
  2. Handle can be a figure handle or a vector of 4 axes handles 
     If Handle=0, no results are plotted
  3. Standard deviation of the BPMs (scalar if all BPMs are the same) 
     These should be in the data file, but this provides an override if not
     found then the default is inf (ie, not used).

  OUTPUTS
  1. For details of the QMS data structure see help quadcenter
     This function added:
     QMS.Offset - Offset computed at each BPM
     QMS.FitParameters - Fit parameter at each BPM
     QMS.FitParametersStd - Sigma of the fit parameter at each BPM
     QMS.BPMStd - BPM sigma at each BPM
     QMS.OffsetSTDMontiCarlo - Monti Carlo estimate of the sigma of the offset (optional)

  2. WarningString = string with warning message if you occurred

0001 function [QMS, WarningString] = quadplotphysics(Input1, FigureHandle, sigmaBPM)
0002 %QUADPLOTPHYSICS - Plots quadrupole centering data in physics units
0003 %  [QMS, WarningString] = quadplot(Input1, Handle, sigmaBPM)
0004 %
0005 %  INPUTS
0006 %  1. Input1 can be a filename for the data or a QMS structure (see help quadcenter for details).
0007 %     If empty or zero inputs, then a dialog box will be provided to select a file.
0008 %  2. Handle can be a figure handle or a vector of 4 axes handles
0009 %     If Handle=0, no results are plotted
0010 %  3. Standard deviation of the BPMs (scalar if all BPMs are the same)
0011 %     These should be in the data file, but this provides an override if not
0012 %     found then the default is inf (ie, not used).
0013 %
0014 %  OUTPUTS
0015 %  1. For details of the QMS data structure see help quadcenter
0016 %     This function added:
0017 %     QMS.Offset - Offset computed at each BPM
0018 %     QMS.FitParameters - Fit parameter at each BPM
0019 %     QMS.FitParametersStd - Sigma of the fit parameter at each BPM
0020 %     QMS.BPMStd - BPM sigma at each BPM
0021 %     QMS.OffsetSTDMontiCarlo - Monti Carlo estimate of the sigma of the offset (optional)
0022 %
0023 %  2. WarningString = string with warning message if you occurred
0024 
0025 %  Written by Greg Portmann
0026 
0027 
0028 % To Do:
0029 % 1. It wouldn't be to difficult to added a LS weight based on slope or even the ideal weighting of std(center).
0030 %    I haven't done it yet because the BPM errors are usually roughly equal at most accelerators.
0031 
0032 
0033 % Remove BPM if it's slope less than MinSlopeFraction * (the maximum slope)
0034 MinSlopeFraction = .25;
0035 
0036 % # of STD of the center calculation allowed
0037 CenterOutlierFactor = 1;
0038 
0039 
0040 QMS = [];
0041 WarningString = '';
0042 
0043 
0044 % Inputs
0045 try
0046     if nargin == 0
0047         [FileName, PathName] = uigetfile('*.mat', 'Select a Quadrupole Center File', [getfamilydata('Directory','DataRoot'), 'QMS', filesep]);
0048         if ~isstr(FileName)
0049             return
0050         else
0051             load([PathName,FileName]);
0052         end
0053     else
0054         if isempty(Input1)
0055             [FileName, PathName] = uigetfile('*.mat', 'Select a Quadrupole Center File', [getfamilydata('Directory','DataRoot'), 'QMS', filesep]);
0056             if ~isstr(FileName)
0057                 return
0058             else
0059                 load([PathName,FileName]);
0060             end
0061         elseif isstr(Input1)
0062             FileName = Input1;
0063             load(FileName);
0064         else
0065             QMS = Input1;
0066             FileName = [];
0067         end
0068     end
0069 catch
0070     error('Problem getting input data');
0071 end
0072 if nargin < 2
0073     FigureHandle = [];
0074 end
0075 
0076 
0077  QMS.QuadraticFit = 0;
0078  
0079 QuadraticFit = QMS.QuadraticFit;       % 0 = linear fit, else quadratic fit
0080 OutlierFactor = QMS.OutlierFactor;     % if abs(data - fit) > OutlierFactor, then remove that BPM
0081 
0082 % Get BPM standard deviation
0083 if nargin < 3
0084     % Get from the data file
0085     if isfield(QMS, 'BPMSTD')
0086         sigmaBPM = QMS.BPMSTD;
0087     else
0088         sigmaBPM = inf;
0089     end
0090 end
0091 if isempty(sigmaBPM)
0092     sigmaBPM = inf;
0093 end
0094 if isnan(sigmaBPM) | isinf(sigmaBPM)
0095     sigmaBPM = inf;
0096     fprintf('   WARNING: BPM standard deviation is unknown, hence there is no BPM outlier condition.\n');
0097 end
0098 sigmaBPM = sigmaBPM(:);
0099 QMS.BPMSTD = sigmaBPM;
0100 
0101 % Get figure handle
0102 if all(FigureHandle ~= 0) 
0103     if isempty(FigureHandle)
0104         FigureHandle = figure;
0105         clf reset
0106         AxesHandles(1) = subplot(3,1,1);
0107         AxesHandles(2) = subplot(3,1,2);
0108         AxesHandles(3) = subplot(3,1,3);
0109         %AxesHandles(4) = subplot(4,1,4);
0110     else
0111         if length(FigureHandle) == 1
0112             FigureHandle = figure(FigureHandle);
0113             clf reset
0114             AxesHandles(1) = subplot(3,1,1);
0115             AxesHandles(2) = subplot(3,1,2);
0116             AxesHandles(3) = subplot(3,1,3);
0117             %AxesHandles(4) = subplot(4,1,4);
0118         elseif length(FigureHandle) == 4
0119             FigureHandle = figure;
0120             clf reset
0121             AxesHandles = FigureHandle;
0122         else
0123             error('Improper size of input FigureHandle');
0124         end
0125     end
0126 end
0127 
0128 Buffer = .01;
0129 HeightBuffer = .08;
0130 if QMS.QuadPlane == 1    
0131     x1 = QMS.x1;
0132     x2 = QMS.x2;
0133     
0134     % Plot setup
0135     if all(FigureHandle ~= 0) 
0136         set(FigureHandle,'units','normal','position',[.0+Buffer .27+Buffer .5-2*Buffer .72-2*Buffer-HeightBuffer]);
0137     end
0138 else
0139     x1 = QMS.y1;
0140     x2 = QMS.y2;
0141     
0142     % Plot setup
0143     if all(FigureHandle ~= 0) 
0144         set(FigureHandle,'units','normal','position',[.5+Buffer .27+Buffer .5-2*Buffer .72-2*Buffer-HeightBuffer]);
0145     end
0146 end
0147 
0148 
0149 [BPMelem1, iNotFound] = findrowindex(QMS.BPMDev, QMS.BPMDevList);
0150 if ~isempty(iNotFound)
0151     error('BPM at the quadrupole not found in the BPM device list');
0152 end
0153 
0154 % Expand sigmaBPM is necessary
0155 if length(sigmaBPM) == 1
0156     sigmaBPM = ones(size(x1,1),1) * sigmaBPM;
0157 end
0158 
0159 N = size(x1,2);
0160 
0161 % Change the number of points
0162 % if 0
0163 %     Ndiv2 = floor(size(x1,2)/2);
0164 %     Npoint1 = Ndiv2;
0165 %     Npoint2 = Ndiv2+2;
0166 %     fprintf('  Using %d points (%d to %d, total %d).', Npoint2-Npoint1+1, Npoint1, Npoint2, N)
0167 %     x1 = x1(:,Npoint1:Npoint2);
0168 %     x2 = x2(:,Npoint1:Npoint2);
0169 %     y1 = y1(:,Npoint1:Npoint2);
0170 %     y2 = y2(:,Npoint1:Npoint2);
0171 %     N = size(x1,2);
0172 %     Ndiv2 = floor(size(x1,2)/2);
0173 % end
0174 
0175 
0176 
0177 
0178 %
0179 % QUAD0 = getquad(QMS, 'Model');
0180 % CM0 = getsp(QMS.CorrFamily, QMS.CorrDevList, 'Model');
0181 %
0182 %
0183 % % Start the corrector a little lower first for hysteresis reasons
0184 % CorrStep = 2 * QMS.CorrDelta / (N-1);
0185 % stepsp(QMS.CorrFamily, -QMS.CorrDelta, QMS.CorrDevList, -1, 'Model');
0186 %
0187 % %XstartModel = getam(BPMxFamily, BPMxDev)
0188 % for i = 1:N
0189 %     % Step the vertical orbit
0190 %     if i ~= 1
0191 %         stepsp(QMS.CorrFamily, CorrStep, QMS.CorrDevList, -1, 'Model');
0192 %     end
0193 %
0194 % %    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);
0195 %
0196 %     %OrbitCorrection(XstartModel, BPMxFamily, BPMxDev, HCMFamily, HCMDev, 2, 'Model');
0197 %
0198 %     if strcmpi(lower(QMS.ModulationMethod), 'sweep')
0199 %         % One dimensional sweep of the quadrupole
0200 %         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0201 %         xm0(:,i) = xm1(:,i);
0202 %         setquad(QMS, i*QMS.QuadDelta+QUAD0, -1, 'Model');
0203 %         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDevList, 'Model');
0204 %     elseif strcmpi(lower(QMS.ModulationMethod), 'bipolar')
0205 %         % Modulate the quadrupole
0206 %         xm0(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0207 %         [xq0(:,i), yq0(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);
0208 %
0209 %         setquad(QMS, QMS.QuadDelta+QUAD0, -1, 'Model');
0210 %         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0211 %         [xq1(:,i), yq1(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);
0212 %
0213 %
0214 %         setquad(QMS,-QMS.QuadDelta+QUAD0, -1, 'Model');
0215 %         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0216 %         [xq2(:,i), yq2(:,i)] = modeltwiss('x', QMS.QuadFamily, QMS.QuadDev);
0217 %
0218 %         setquad(QMS, QUAD0, -1, 'Model');
0219 %     elseif strcmpi(lower(QMS.ModulationMethod), 'unipolar')
0220 %         % Modulate the quadrupole
0221 %         xm1(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0222 %         xm0(:,i) = x1(:,i);
0223 %         setquad(QMS, QMS.QuadDelta+QUAD0, -1, 'Model');
0224 %         xm2(:,i) = getam(QMS.BPMFamily, QMS.BPMDev, 'Model');
0225 %         setquad(QMS, QUAD0, -1, 'Model');
0226 %     end
0227 % end
0228 %
0229 % setquad(QMS, QUAD0, -1, 'Model');
0230 % setsp(QMS.CorrFamily, CM0, QMS.CorrDevList, -1, 'Model');
0231 %
0232 % xq0 = 1000*xq0;
0233 % xq1 = 1000*xq1;
0234 % xq2 = 1000*xq2;
0235 %
0236 % yq0 = 1000*yq0;
0237 % yq1 = 1000*yq1;
0238 % yq2 = 1000*yq2;
0239 %
0240 % x = xm0 + yq0-(xm0-xm0(3));
0241 % x = xm0(3) + yq0-yq0(3);
0242 % x = x';
0243 
0244 
0245 % Convert to physics units
0246 %x0 = hw2physics('BPMx', 'Monitor', x0, QMS.BPMDevList);
0247 %x1 = hw2physics('BPMx', 'Monitor', x1, QMS.BPMDevList);
0248 %x2 = hw2physics('BPMx', 'Monitor', x2, QMS.BPMDevList);
0249 
0250 
0251 Gx = getgain('BPMx', QMS.BPMDevList);
0252 Gy = getgain('BPMy', QMS.BPMDevList);
0253 C = getcrunch('BPMx', QMS.BPMDevList);
0254 R = getroll('BPMx', QMS.BPMDevList);
0255 
0256 for i = 1:length(Gx)
0257     m = [cos(R(i)) -sin(R(i)); sin(R(i)) cos(R(i))] * [1 C(i); C(i) 1] * [Gx(i) 0;0 Gy(i)] / sqrt(1-C(i)^2);
0258 
0259     for j = 1:size(QMS.x1,2)
0260         if j == 1
0261             if isfield(QMS, 'x0')
0262                 a = m * [QMS.x0(i,j); QMS.y0(i,j)];
0263                 QMS.x0(i,j) = a(1);
0264                 QMS.y0(i,j) = a(2);
0265             end
0266         end
0267 
0268         a = m * [QMS.x1(i,j); QMS.y1(i,j)];
0269         QMS.x1(i,j) = a(1);
0270         QMS.y1(i,j) = a(2);
0271         
0272         a = m * [QMS.x2(i,j); QMS.y2(i,j)];
0273         QMS.x2(i,j) = a(1);
0274         QMS.y2(i,j) = a(2);
0275     end
0276 end
0277 
0278 x1 = QMS.y1;
0279 x2 = QMS.y2; 
0280 
0281 
0282 
0283 if isfield(QMS, 'OldCenter')
0284  %   QMS.OldCenter = hw2physics(QMS.BPMFamily, 'Monitor', QMS.OldCenter, QMS.BPMDev);
0285 end
0286 
0287 if isfield(QMS, 'Orbit0')
0288   %  QMS.Orbit0 = hw2physics(QMS.Orbit0);
0289     BPMUnitsString = QMS.Orbit0.UnitsString;
0290 else
0291     BPMUnitsString = 'mm';
0292 end
0293 
0294 
0295 % Fit verses the position at the BPM next to the quadrupole
0296 if strcmpi(QMS.ModulationMethod, 'sweep')
0297     % One dimensional sweep of the quadrupole
0298     %x = x1(BPMelem1,:)';
0299     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
0300 elseif strcmpi(QMS.ModulationMethod, 'bipolar')
0301     % Modulation of the quadrupole
0302     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
0303 elseif strcmpi(QMS.ModulationMethod, 'unipolar')
0304     % Unipolar modulation of the quadrupole
0305     %x = x1(BPMelem1,:)';
0306     x = (x1(BPMelem1,:)' + x2(BPMelem1,:)')/2;     
0307 end
0308 
0309 
0310 % Figure #1
0311 if all(FigureHandle ~= 0) 
0312     axes(AxesHandles(1));
0313     %plot(linspace(-DelHCM,DelHCM,3), x2-x1);
0314     plot(x, x2-x1, '-x');
0315     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0316     ylabel(sprintf('%s [%s]', QMS.BPMFamily, BPMUnitsString));
0317     %plot(1000*x, 1000*(x2-x1), '-x');
0318     %xlabel(sprintf('%s(%d,%d), raw values [mm]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2)));
0319     %ylabel(sprintf('%s [mm]', QMS.BPMFamily));
0320     if isempty(FileName)
0321         title(sprintf('Center for %s(%d,%d) %s(%d,%d)', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev), 'interpreter', 'none');
0322     else
0323         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');
0324     end
0325     grid on
0326     axis tight
0327 end
0328 
0329 if isempty(FileName)
0330     fprintf('   Calculating the center of %s(%d,%d) using %s(%d,%d)\n', QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev);
0331 else
0332     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);
0333 end
0334 fprintf('   Quadrupole modulation delta = %.3f amps, max. corrector step = %.3f amps\n', QMS.QuadDelta, QMS.CorrDelta);
0335 
0336 
0337 % Least squares fit
0338 merit = x2 - x1;
0339 if QuadraticFit
0340     X = [ones(size(x)) x x.^2];
0341     fprintf('   %d point parabolic least squares fit\n', N);
0342 else
0343     X = [ones(size(x)) x];
0344     fprintf('   %d point linear least squares fit\n', N);
0345 end
0346 
0347 
0348 % Axes #2
0349 if all(FigureHandle ~= 0) 
0350     axes(AxesHandles(2));
0351     xx = linspace(x(1), x(end), 200);
0352 end
0353 
0354 iBPMOutlier = [];
0355 invXX   = inv(X'*X);
0356 invXX_X = invXX * X';
0357 
0358 for i = 1:size(x1,1)
0359     % least-square fit: m = slope and b = Y-intercept
0360     b = invXX_X * merit(i,:)';
0361     bhat(i,:) = b';
0362     
0363     % Should equal
0364     %b = X \merit(i,:)';
0365     %bhat1(i,:) = b';
0366     
0367     % Standard deviation
0368     bstd = sigmaBPM(i) * invXX; 
0369     bhatstd(i,:) = diag(bstd)';  % hopefully cross-correlation terms are small
0370     
0371     if all(FigureHandle ~= 0)   
0372         if QuadraticFit
0373             y = b(3)*xx.^2 + b(2)*xx + b(1);
0374         else
0375             y = b(2)*xx + b(1);
0376         end
0377 %        plot(xx, y); hold on
0378     end
0379     
0380     % Outlier condition: remove if the error between the fit and the data is greater than OutlierFactor
0381     if QuadraticFit
0382         y = b(3)*x.^2 + b(2)*x + b(1);
0383     else
0384         y = b(2)*x + b(1);
0385     end
0386     if max(abs(y - merit(i,:)')) > OutlierFactor * sigmaBPM(i)    % OutlierFactor was absolute max(abs(y - merit(i,:)')) > OutlierFactor
0387         iBPMOutlier = [iBPMOutlier;i];
0388     end
0389     
0390     if QuadraticFit
0391         % Quadratic fit
0392         if b(2) > 0
0393             offset(i,1) = (-b(2) + sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
0394         else
0395             offset(i,1) = (-b(2) - sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
0396         end
0397         if ~isreal(offset(i,1))
0398             % (b^2-4ac) can be negative but it will only happen if the slope is very small.  The offset
0399             % should just get thrown out later as an outlier but change the solution to the minimum of the parabola.
0400             offset(i,1) = -b(2) / b(1) / 2;
0401         end
0402     else
0403         % Linear fit
0404         offset(i,1) = -b(1)/b(2); 
0405     end
0406 end
0407 
0408 QMS.Offset = offset;
0409 QMS.FitParameters = bhat;
0410 QMS.FitParametersStd = bhatstd;
0411 QMS.BPMStd = sigmaBPM;
0412 
0413 
0414 % % Label axes #2
0415 % if all(FigureHandle ~= 0)
0416 %     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0417 %     ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString));
0418 %     grid on
0419 %     axis tight
0420 % end
0421 
0422 
0423 % Compute offset for different conditions
0424 fprintf('   %d total BPMs\n', length(offset));
0425 fprintf('   BPMs are removed for 1. Bad Status, 2. BPM Outlier, 3. Small Slope, or 4. Center Outlier\n');
0426 %m1 = mean(offset);
0427 %s1 = std(offset);
0428 %fprintf('   0. Mean = %.5f %s, STD = %.5f %s, all %d BPMs\n', m1, BPMUnitsString, s1, BPMUnitsString, length(offset));
0429 
0430 
0431 % Remove bad Status BPMs
0432 iStatus = find(QMS.BPMStatus==0);
0433 iBad = iStatus;
0434 if length(iBad) == length(offset)
0435     error('All the BPMs have bad status');
0436 end
0437 offset1 = offset;
0438 offset1(iBad) = [];
0439 m2 = mean(offset1);
0440 s2 = std(offset1);
0441 fprintf('   1. Mean = %+.5f %s, STD = %.5f %s, %2d points with bad status\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad));
0442 
0443 % Remove bad Status + Outliers
0444 iBad = unique([iStatus; iBPMOutlier]);
0445 if length(iBad) == length(offset)
0446     error('All BPMs either have bad status or failed the BPM outlier condition');
0447 end
0448 offset1 = offset;
0449 offset1(iBad) = [];
0450 m2 = mean(offset1);
0451 s2 = std(offset1);
0452 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);
0453 
0454 
0455 % Remove bad Status + Small slopes
0456 %iSlope = find(abs(bhat(:,2)) < max(abs(bhat(:,2)))*MinSlopeFraction);
0457 
0458 % Look for slope outliers
0459 Slopes = abs(bhat(:,2));
0460 [Slopes, i] = sort(Slopes);
0461 Slopes = Slopes(round(end/2):end);  % remove the first half
0462 if length(Slopes) > 5
0463     SlopesMax = Slopes(end-4); 
0464 else
0465     SlopesMax = Slopes(end); 
0466 end
0467 %i = find(abs(Slopes-mean(Slopes)) > 2 * std(Slopes));
0468 %Slopes(i) = [];
0469 
0470 iSlope = find(abs(bhat(:,2)) < SlopesMax * MinSlopeFraction);
0471 iBad = unique([iStatus; iBPMOutlier; iSlope]);
0472 if length(iBad) == length(offset)
0473     error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition');
0474 end
0475 offset1 = offset;
0476 offset1(iBad) = [];
0477 m2 = mean(offset1);
0478 s2 = std(offset1);
0479 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);
0480 
0481 
0482 % Offset outlier offsets-mean(offsets) greater than 1 std
0483 itotal = (1:length(offset))';
0484 iok = itotal;
0485 
0486 offset1 = offset;
0487 offset1(iBad) = [];
0488 iok(iBad) = [];
0489 
0490 i = find(abs(offset1-mean(offset1)) > CenterOutlierFactor * std(offset1));
0491 iCenterOutlier = iok(i);
0492 iBad = unique([iBad; iCenterOutlier]);
0493 if length(iBad) == length(offset)
0494     error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition, , or failed the center outlier condition');
0495 end
0496 offset1(i) = [];
0497 iok(i) = [];
0498 
0499 m2 = mean(offset1);
0500 s2 = std(offset1);
0501 QMS.GoodIndex = iok;
0502 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);
0503 
0504 
0505 NN = length(offset);
0506 
0507 % Axes #4
0508 if all(FigureHandle ~= 0) 
0509     axes(AxesHandles(3));
0510     [xx1,yy1]=stairs(1:NN,offset);
0511     offset1 = offset;
0512     offset1(iBad) = NaN*ones(length(iBad),1);
0513     [xx2, yy2] = stairs(1:NN, offset1);
0514     plot(xx1,yy1,'r', xx2,yy2,'b');
0515     ylabel(sprintf('BPM Center [%s]', BPMUnitsString));
0516     %plot(xx1,1000*yy1,'r', xx2,1000*yy2,'b');
0517     %ylabel(sprintf('BPM Center [mm]'));
0518     xlabel('BPM Number');
0519     grid on
0520     axis tight
0521     %xaxis([1 NN+1]);
0522     %axis([1 NN+1 min(offset1)-.1e-6 max(offset1)+.1e-6]);
0523     %axis([1 NN+1 1000*[min(offset1)-.1e-6 max(offset1)+.1e-6]]);
0524     axis([1 NN+1 min(offset1)-.1 max(offset1)+.1]);
0525 end
0526 
0527 
0528 % Axes #2
0529 
0530 if all(FigureHandle ~= 0) 
0531     if 0
0532         % Plot red line over the bad lines
0533         axes(AxesHandles(2));
0534         for j = 1:length(iBad)
0535             i = iBad(j);
0536             if QuadraticFit
0537                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0538             else
0539                 y = bhat(i,2)*xx + bhat(i,1);
0540             end
0541             plot(xx, y,'r'); 
0542         end
0543         hold off
0544         axis tight
0545     else
0546         % Only plot the good data
0547         axes(AxesHandles(2));
0548         yy = [];
0549         for i = 1:size(x1,1)
0550             if ~any(i == iBad)
0551                 if QuadraticFit
0552                     y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0553                 else
0554                     y = bhat(i,2)*xx + bhat(i,1);
0555                 end
0556                 yy = [yy;y];
0557             end
0558         end
0559         %plot(xx, yy);
0560         plot(1000*xx, 1000*yy); 
0561         hold off
0562         grid on
0563         axis tight
0564     end
0565     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0566     ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString));
0567     %xlabel(sprintf('%s(%d,%d), raw values [mm]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2)));
0568     %ylabel(sprintf('BPM LS Fit [mm]'));
0569     grid on
0570     axis tight
0571 end
0572 
0573 
0574 
0575 if ~isempty(iStatus)
0576     if ~isempty(find(iStatus==BPMelem1))
0577         fprintf('   WARNING: BPM(%d,%d) has a bad status\n', QMS.BPMDev(1), QMS.BPMDev(2));
0578         WarningString = sprintf('BPM(%d,%d) has a bad status', QMS.BPMDev(1), QMS.BPMDev(2));
0579     end
0580 end
0581 if ~isempty(iBPMOutlier)
0582     if ~isempty(find(iBPMOutlier==BPMelem1))
0583         fprintf('   WARNING: BPM(%d,%d) removed due to outlier (based on std(BPM))\n', QMS.BPMDev(1), QMS.BPMDev(2));
0584         WarningString = sprintf('BPM(%d,%d) removed due to outlier (based on std(BPM))', QMS.BPMDev(1), QMS.BPMDev(2));
0585     end
0586 end
0587 if ~isempty(iSlope)
0588     if ~isempty(find(iSlope==BPMelem1))
0589         fprintf('   WARNING: BPM(%d,%d) slope is too small\n', QMS.BPMDev(1), QMS.BPMDev(2));
0590         WarningString = sprintf('BPM(%d,%d) slope is too small', QMS.BPMDev(1), QMS.BPMDev(2));
0591     end
0592 end
0593 if ~isempty(iCenterOutlier)
0594     if ~isempty(find(iCenterOutlier==BPMelem1))
0595         fprintf('   WARNING: BPM(%d,%d) removed due to outlier (based on all the centers)\n', QMS.BPMDev(1), QMS.BPMDev(2));
0596         WarningString = sprintf('BPM(%d,%d) ', QMS.BPMDev(1), QMS.BPMDev(2));
0597     end
0598 end
0599 
0600 
0601 % % Axes #3
0602 % if all(FigureHandle ~= 0)
0603 %     axes(AxesHandles(3));
0604 %     iii=1:NN;
0605 %     iii(iBad)=[];
0606 %     for j = 1:length(iii)
0607 %         i = iii(j);
0608 %
0609 %         if 1
0610 %             % Plot fit
0611 %             if QuadraticFit
0612 %                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0613 %             else
0614 %                 y = bhat(i,2)*xx + bhat(i,1);
0615 %             end
0616 %             if all(FigureHandle ~= 0)
0617 %                 plot(xx, y,'b'); hold on
0618 %             end
0619 %         else
0620 %             % Plot error in fit
0621 %             if QuadraticFit
0622 %                 y = bhat(i,3)*x.^2 + bhat(i,2)*x + bhat(i,1);
0623 %             else
0624 %                 y = bhat(i,2)*x + bhat(i,1);
0625 %             end
0626 %             plot(x, y - merit(i,:)','b'); hold on
0627 %         end
0628 %     end
0629 %     hold off
0630 %     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0631 %     ylabel(sprintf('Final %s Merit Fun [%s]', QMS.BPMFamily, BPMUnitsString));
0632 %     grid on
0633 %     axis tight
0634 %     orient tall
0635 % end
0636 
0637 if ~isempty(QMS.OldCenter)
0638     fprintf('   Starting Offset %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, QMS.OldCenter, BPMUnitsString);
0639 end
0640 fprintf('   New Offset      %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, m2, BPMUnitsString);
0641 
0642 QMS.Center = m2;
0643 QMS.CenterSTD = s2;
0644 
0645 
0646 if all(FigureHandle ~= 0)
0647     addlabel(1, 0, datestr(QMS.TimeStamp));
0648     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)));
0649 end
0650 
0651 fprintf('\n');
0652 
0653 
0654 % Get and plot errors
0655 if all(FigureHandle ~= 0)
0656     QMS = quaderrors(QMS, gcf+1);
0657 else
0658     % This is a cluge for now
0659     QMS = quaderrors(QMS, 0);
0660 end
0661 
0662 
0663 %QMS = orderfields(QMS);