0001 function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag)
0002
0003
0004
0005
0006
0007
0008 Buffer = .03;
0009 HeightBuffer = .08;
0010
0011
0012
0013 MontiCarloFlag = 1;
0014 [ST,I] = dbstack;
0015 for i = 1:length(ST)
0016 if strcmpi(ST(i).file,'quadplot.m')
0017 MontiCarloFlag = 0;
0018 end
0019 end
0020
0021
0022
0023 try
0024 WarningString = [];
0025 if nargin == 0
0026 [FileName, PathName] = uigetfile('*.mat', 'Input Quadrupole Center file.');
0027 if ~isstr(FileName)
0028 return
0029 else
0030 eval(['load ',PathName,FileName]);
0031 end
0032 else
0033 if isstr(Input1)
0034 FileName = Input1;
0035 eval(['load ', FileName]);
0036 elseif isempty(Input1)
0037 [FileName, PathName] = uigetfile('*.mat', 'Input Quadrupole Center file.');
0038 if ~isstr(FileName)
0039 return
0040 else
0041 eval(['load ',PathName,FileName]);
0042 end
0043 else
0044 QMS = Input1;
0045 FileName = [];
0046 end
0047 end
0048 catch
0049 error('Problem getting input data');
0050 end
0051
0052
0053
0054 if nargin < 2
0055 FigureHandle = gcf;
0056 else
0057 if all(FigureHandle ~= 0)
0058 figure(FigureHandle);
0059 end
0060 end
0061
0062 if all(FigureHandle ~= 0)
0063 clf reset
0064 if MontiCarloFlag
0065 AxesHandles(1) = subplot(3,1,1);
0066 AxesHandles(2) = subplot(3,1,2);
0067 AxesHandles(3) = subplot(3,1,3);
0068 else
0069 AxesHandles(1) = subplot(2,1,1);
0070 AxesHandles(2) = subplot(2,1,2);
0071 end
0072
0073 if QMS.QuadPlane ~= 1
0074 set(FigureHandle,'units','normal','position',[Buffer .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]);
0075
0076 else
0077 set(FigureHandle,'units','normal','position',[.5+.003 .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]);
0078
0079 end
0080
0081
0082 axes(AxesHandles(1));
0083 errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),'b');
0084 hold on;
0085 errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),'.r');
0086 errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,1), QMS.FitParametersStd(QMS.GoodIndex,1),'.b');
0087 hold off
0088 title('Linear Fit');
0089 ylabel('Y-Intercept');
0090
0091 grid
0092 xaxis([0 length(QMS.FitParameters(:,1))+1]);
0093
0094
0095
0096 axes(AxesHandles(2));
0097 errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),'b');
0098 hold on
0099 errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),'.r');
0100 errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,2), QMS.FitParametersStd(QMS.GoodIndex,2),'.b');
0101 hold off
0102 xlabel('BPM Number');
0103 ylabel('Slope');
0104
0105 grid
0106 xaxis([0 length(QMS.FitParameters(:,1))+1]);
0107 end
0108
0109
0110 if MontiCarloFlag
0111
0112 N = 20000;
0113 QMS.OffsetMean = NaN*ones(length(QMS.FitParameters(:,2)),1);
0114 QMS.OffsetStd = NaN*ones(length(QMS.FitParameters(:,2)),1);
0115
0116 OffsetMean1 = QMS.OffsetMean;
0117 OffsetStd1 = QMS.OffsetStd;
0118 NormalRV1 = randn(N,1);
0119 NormalRV2 = randn(N,1);
0120 if QMS.QuadraticFit == 1
0121 NormalRV3 = randn(N,1);
0122 end
0123 for i = QMS.GoodIndex(:)'
0124 rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1;
0125 rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2;
0126 if QMS.QuadraticFit == 1
0127 rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3;
0128 x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic);
0129 else
0130 x = -rvinter ./ rvslope;
0131 end
0132 OffsetMean1(i) = mean(x);
0133 OffsetStd1(i) = std(x);
0134 end
0135
0136 OffsetMean2 = QMS.OffsetMean;
0137 OffsetStd2 = QMS.OffsetStd;
0138 NormalRV1 = randn(N,1);
0139 NormalRV2 = randn(N,1);
0140 if QMS.QuadraticFit == 1
0141 NormalRV3 = randn(N,1);
0142 end
0143 for i = QMS.GoodIndex(:)'
0144 rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1;
0145 rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2;
0146 if QMS.QuadraticFit == 1
0147 rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3;
0148 x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic);
0149 else
0150 x = -rvinter ./ rvslope;
0151 end
0152 OffsetMean2(i) = mean(x);
0153 OffsetStd2(i) = std(x);
0154 end
0155
0156 QMS.OffsetMean(:,1) = [OffsetMean1 + OffsetMean2] / 2;
0157 QMS.OffsetStd(:,1) = [OffsetStd1 + OffsetStd2] / 2;
0158
0159
0160 if all(FigureHandle ~= 0)
0161
0162
0163
0164
0165 axes(AxesHandles(3));
0166
0167
0168 errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), QMS.OffsetStd(QMS.GoodIndex),'.b');
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183 hold off
0184 ylabel('X-Intercept = - (Y-Intercept)/Slope');
0185 ylabel('Monte Carlo');
0186 xlabel('BPM Number');
0187 grid
0188 xaxis([0 length(QMS.Offset)+1]);
0189
0190 end
0191
0192 QMS.OffsetSTDMontiCarlo = sqrt((sum(QMS.OffsetStd(QMS.GoodIndex).^2))/length(QMS.GoodIndex));
0193 end