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

  Written by Greg Portmann

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 = 'm';
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(x, x2-x1, '-x');
0314     %plot(linspace(-DelHCM,DelHCM,3), x2-x1);
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     if isempty(FileName)
0318         title(sprintf('Center for %s(%d,%d) %s(%d,%d)', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), QMS.QuadFamily, QMS.QuadDev), 'interpreter', 'none');
0319     else
0320         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');
0321     end
0322     grid on
0323     axis tight
0324 end
0325 
0326 if isempty(FileName)
0327     fprintf('   Calculating the center of %s(%d,%d) using %s(%d,%d)\n', QMS.QuadFamily, QMS.QuadDev, QMS.BPMFamily, QMS.BPMDev);
0328 else
0329     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);
0330 end
0331 fprintf('   Quadrupole modulation delta = %.3f amps, max. corrector step = %.3f amps\n', QMS.QuadDelta, QMS.CorrDelta);
0332 
0333 
0334 % Least squares fit
0335 merit = x2 - x1;
0336 if QuadraticFit
0337     X = [ones(size(x)) x x.^2];
0338     fprintf('   %d point parabolic least squares fit\n', N);
0339 else
0340     X = [ones(size(x)) x];
0341     fprintf('   %d point linear least squares fit\n', N);
0342 end
0343 
0344 
0345 % Axes #2
0346 if all(FigureHandle ~= 0) 
0347     axes(AxesHandles(2));
0348     xx = linspace(x(1), x(end), 200);
0349 end
0350 
0351 iBPMOutlier = [];
0352 invXX   = inv(X'*X);
0353 invXX_X = invXX * X';
0354 
0355 for i = 1:size(x1,1)
0356     % least-square fit: m = slope and b = Y-intercept
0357     b = invXX_X * merit(i,:)';
0358     bhat(i,:) = b';
0359     
0360     % Should equal
0361     %b = X \merit(i,:)';
0362     %bhat1(i,:) = b';
0363     
0364     % Standard deviation
0365     bstd = sigmaBPM(i) * invXX; 
0366     bhatstd(i,:) = diag(bstd)';  % hopefully cross-correlation terms are small
0367     
0368     if all(FigureHandle ~= 0)   
0369         if QuadraticFit
0370             y = b(3)*xx.^2 + b(2)*xx + b(1);
0371         else
0372             y = b(2)*xx + b(1);
0373         end
0374 %        plot(xx, y); hold on
0375     end
0376     
0377     % Outlier condition: remove if the error between the fit and the data is greater than OutlierFactor
0378     if QuadraticFit
0379         y = b(3)*x.^2 + b(2)*x + b(1);
0380     else
0381         y = b(2)*x + b(1);
0382     end
0383     if max(abs(y - merit(i,:)')) > OutlierFactor * sigmaBPM(i)    % OutlierFactor was absolute max(abs(y - merit(i,:)')) > OutlierFactor
0384         iBPMOutlier = [iBPMOutlier;i];
0385     end
0386     
0387     if QuadraticFit
0388         % Quadratic fit
0389         if b(2) > 0
0390             offset(i,1) = (-b(2) + sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
0391         else
0392             offset(i,1) = (-b(2) - sqrt(b(2)^2 - 4*b(3)*b(1))) / (2*b(3));
0393         end
0394         if ~isreal(offset(i,1))
0395             % (b^2-4ac) can be negative but it will only happen if the slope is very small.  The offset
0396             % should just get thrown out later as an outlier but change the solution to the minimum of the parabola.
0397             offset(i,1) = -b(2) / b(1) / 2;
0398         end
0399     else
0400         % Linear fit
0401         offset(i,1) = -b(1)/b(2); 
0402     end
0403 end
0404 
0405 QMS.Offset = offset;
0406 QMS.FitParameters = bhat;
0407 QMS.FitParametersStd = bhatstd;
0408 QMS.BPMStd = sigmaBPM;
0409 
0410 
0411 % % Label axes #2
0412 % if all(FigureHandle ~= 0)
0413 %     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0414 %     ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString));
0415 %     grid on
0416 %     axis tight
0417 % end
0418 
0419 
0420 % Compute offset for different conditions
0421 fprintf('   %d total BPMs\n', length(offset));
0422 fprintf('   BPMs are removed for 1. Bad Status, 2. BPM Outlier, 3. Small Slope, or 4. Center Outlier\n');
0423 %m1 = mean(offset);
0424 %s1 = std(offset);
0425 %fprintf('   0. Mean = %.5f %s, STD = %.5f %s, all %d BPMs\n', m1, BPMUnitsString, s1, BPMUnitsString, length(offset));
0426 
0427 
0428 % Remove bad Status BPMs
0429 iStatus = find(QMS.BPMStatus==0);
0430 iBad = iStatus;
0431 if length(iBad) == length(offset)
0432     error('All the BPMs have bad status');
0433 end
0434 offset1 = offset;
0435 offset1(iBad) = [];
0436 m2 = mean(offset1);
0437 s2 = std(offset1);
0438 fprintf('   1. Mean = %+.5f %s, STD = %.5f %s, %2d points with bad status\n', m2, BPMUnitsString, s2, BPMUnitsString, length(iBad));
0439 
0440 % Remove bad Status + Outliers
0441 iBad = unique([iStatus; iBPMOutlier]);
0442 if length(iBad) == length(offset)
0443     error('All BPMs either have bad status or failed the BPM outlier condition');
0444 end
0445 offset1 = offset;
0446 offset1(iBad) = [];
0447 m2 = mean(offset1);
0448 s2 = std(offset1);
0449 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);
0450 
0451 
0452 % Remove bad Status + Small slopes
0453 %iSlope = find(abs(bhat(:,2)) < max(abs(bhat(:,2)))*MinSlopeFraction);
0454 
0455 % Look for slope outliers
0456 Slopes = abs(bhat(:,2));
0457 [Slopes, i] = sort(Slopes);
0458 Slopes = Slopes(round(end/2):end);  % remove the first half
0459 if length(Slopes) > 5
0460     SlopesMax = Slopes(end-4); 
0461 else
0462     SlopesMax = Slopes(end); 
0463 end
0464 %i = find(abs(Slopes-mean(Slopes)) > 2 * std(Slopes));
0465 %Slopes(i) = [];
0466 
0467 iSlope = find(abs(bhat(:,2)) < SlopesMax * MinSlopeFraction);
0468 iBad = unique([iStatus; iBPMOutlier; iSlope]);
0469 if length(iBad) == length(offset)
0470     error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition');
0471 end
0472 offset1 = offset;
0473 offset1(iBad) = [];
0474 m2 = mean(offset1);
0475 s2 = std(offset1);
0476 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);
0477 
0478 
0479 % Offset outlier offsets-mean(offsets) greater than 1 std
0480 itotal = (1:length(offset))';
0481 iok = itotal;
0482 
0483 offset1 = offset;
0484 offset1(iBad) = [];
0485 iok(iBad) = [];
0486 
0487 i = find(abs(offset1-mean(offset1)) > CenterOutlierFactor * std(offset1));
0488 iCenterOutlier = iok(i);
0489 iBad = unique([iBad; iCenterOutlier]);
0490 if length(iBad) == length(offset)
0491     error('All BPMs either have bad status, failed the BPM outlier condition, or failed the slope condition, , or failed the center outlier condition');
0492 end
0493 offset1(i) = [];
0494 iok(i) = [];
0495 
0496 m2 = mean(offset1);
0497 s2 = std(offset1);
0498 QMS.GoodIndex = iok;
0499 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);
0500 
0501 
0502 NN = length(offset);
0503 
0504 % Axes #4
0505 if all(FigureHandle ~= 0) 
0506     axes(AxesHandles(3));
0507     [xx1,yy1]=stairs(1:NN,offset);
0508     offset1 = offset;
0509     offset1(iBad) = NaN*ones(length(iBad),1);
0510     [xx2, yy2] = stairs(1:NN, offset1);
0511     plot(xx1,yy1,'r', xx2,yy2,'b');
0512     xlabel('BPM Number');
0513     ylabel(sprintf('BPM Center [%s]', BPMUnitsString));
0514     grid on
0515     axis tight
0516     %xaxis([0 NN+1]);
0517     axis([0 NN+1 min(offset1)-.1e-6 max(offset1)+.1e-6]);
0518 end
0519 
0520 
0521 % Axes #2
0522 
0523 if all(FigureHandle ~= 0) 
0524     if 0
0525         % Plot red line over the bad lines
0526         axes(AxesHandles(2));
0527         for j = 1:length(iBad)
0528             i = iBad(j);
0529             if QuadraticFit
0530                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0531             else
0532                 y = bhat(i,2)*xx + bhat(i,1);
0533             end
0534             plot(xx, y,'r'); 
0535         end
0536         hold off
0537         axis tight
0538     else
0539         % Only plot the good data
0540         axes(AxesHandles(2));
0541         yy = [];
0542         for i = 1:size(x1,1)
0543             if ~any(i == iBad)
0544                 if QuadraticFit
0545                     y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0546                 else
0547                     y = bhat(i,2)*xx + bhat(i,1);
0548                 end
0549                 yy = [yy;y];
0550             end
0551         end
0552         plot(xx, yy); 
0553         hold off
0554         grid on
0555         axis tight
0556     end
0557     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0558     ylabel(sprintf('BPM LS Fit [%s]', BPMUnitsString));
0559     grid on
0560     axis tight
0561 end
0562 
0563 
0564 
0565 if ~isempty(iStatus)
0566     if ~isempty(find(iStatus==BPMelem1))
0567         fprintf('   WARNING: BPM(%d,%d) has a bad status\n', QMS.BPMDev(1), QMS.BPMDev(2));
0568         WarningString = sprintf('BPM(%d,%d) has a bad status', QMS.BPMDev(1), QMS.BPMDev(2));
0569     end
0570 end
0571 if ~isempty(iBPMOutlier)
0572     if ~isempty(find(iBPMOutlier==BPMelem1))
0573         fprintf('   WARNING: BPM(%d,%d) removed due to outlier (based on std(BPM))\n', QMS.BPMDev(1), QMS.BPMDev(2));
0574         WarningString = sprintf('BPM(%d,%d) removed due to outlier (based on std(BPM))', QMS.BPMDev(1), QMS.BPMDev(2));
0575     end
0576 end
0577 if ~isempty(iSlope)
0578     if ~isempty(find(iSlope==BPMelem1))
0579         fprintf('   WARNING: BPM(%d,%d) slope is too small\n', QMS.BPMDev(1), QMS.BPMDev(2));
0580         WarningString = sprintf('BPM(%d,%d) slope is too small', QMS.BPMDev(1), QMS.BPMDev(2));
0581     end
0582 end
0583 if ~isempty(iCenterOutlier)
0584     if ~isempty(find(iCenterOutlier==BPMelem1))
0585         fprintf('   WARNING: BPM(%d,%d) removed due to outlier (based on all the centers)\n', QMS.BPMDev(1), QMS.BPMDev(2));
0586         WarningString = sprintf('BPM(%d,%d) ', QMS.BPMDev(1), QMS.BPMDev(2));
0587     end
0588 end
0589 
0590 
0591 % % Axes #3
0592 % if all(FigureHandle ~= 0)
0593 %     axes(AxesHandles(3));
0594 %     iii=1:NN;
0595 %     iii(iBad)=[];
0596 %     for j = 1:length(iii)
0597 %         i = iii(j);
0598 %
0599 %         if 1
0600 %             % Plot fit
0601 %             if QuadraticFit
0602 %                 y = bhat(i,3)*xx.^2 + bhat(i,2)*xx + bhat(i,1);
0603 %             else
0604 %                 y = bhat(i,2)*xx + bhat(i,1);
0605 %             end
0606 %             if all(FigureHandle ~= 0)
0607 %                 plot(xx, y,'b'); hold on
0608 %             end
0609 %         else
0610 %             % Plot error in fit
0611 %             if QuadraticFit
0612 %                 y = bhat(i,3)*x.^2 + bhat(i,2)*x + bhat(i,1);
0613 %             else
0614 %                 y = bhat(i,2)*x + bhat(i,1);
0615 %             end
0616 %             plot(x, y - merit(i,:)','b'); hold on
0617 %         end
0618 %     end
0619 %     hold off
0620 %     xlabel(sprintf('%s(%d,%d), raw values [%s]', QMS.BPMFamily, QMS.BPMDev(1), QMS.BPMDev(2), BPMUnitsString));
0621 %     ylabel(sprintf('Final %s Merit Fun [%s]', QMS.BPMFamily, BPMUnitsString));
0622 %     grid on
0623 %     axis tight
0624 %     orient tall
0625 % end
0626 
0627 if ~isempty(QMS.OldCenter)
0628     fprintf('   Starting Offset %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, QMS.OldCenter, BPMUnitsString);
0629 end
0630 fprintf('   New Offset      %s(%d,%d) = %+f [%s]\n', QMS.BPMFamily, QMS.BPMDev, m2, BPMUnitsString);
0631 
0632 QMS.Center = m2;
0633 QMS.CenterSTD = s2;
0634 
0635 
0636 if all(FigureHandle ~= 0)
0637     addlabel(1, 0, datestr(QMS.TimeStamp));
0638     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)));
0639 end
0640 
0641 fprintf('\n');
0642 
0643 
0644 % Get and plot errors
0645 if all(FigureHandle ~= 0)
0646     QMS = quaderrors(QMS, gcf+1);
0647 else
0648     % This is a cluge for now
0649     QMS = quaderrors(QMS, 0);
0650 end
0651 
0652 
0653 %QMS = orderfields(QMS);