1 | function QMS = quaderrors(Input1, FigureHandle, MontiCarloFlag) |
---|
2 | %QUADERRORS - Computes the error bars for quadrupole center measurement |
---|
3 | % QMS = quaderrors(FileName, FigureHandle, MontiCarloFlag) |
---|
4 | % |
---|
5 | % See Also quadplot, quadcenter |
---|
6 | |
---|
7 | % |
---|
8 | %% Written by Greg Portmann |
---|
9 | |
---|
10 | |
---|
11 | Buffer = .03; |
---|
12 | HeightBuffer = .08; |
---|
13 | |
---|
14 | |
---|
15 | % Don't do the monticarlo plot if being called from quadplot |
---|
16 | MontiCarloFlag = 1; |
---|
17 | [ST,I] = dbstack; |
---|
18 | for i = 1:length(ST) |
---|
19 | if strcmpi(ST(i).file,'quadplot.m') |
---|
20 | MontiCarloFlag = 0; |
---|
21 | end |
---|
22 | end |
---|
23 | |
---|
24 | |
---|
25 | % Inputs |
---|
26 | try |
---|
27 | WarningString = []; |
---|
28 | if nargin == 0 |
---|
29 | [FileName, PathName] = uigetfile('*.mat', 'Input Quadrupole Center file.'); |
---|
30 | if ~isstr(FileName) |
---|
31 | return |
---|
32 | else |
---|
33 | eval(['load ',PathName,FileName]); |
---|
34 | end |
---|
35 | else |
---|
36 | if isstr(Input1) |
---|
37 | FileName = Input1; |
---|
38 | eval(['load ', FileName]); |
---|
39 | elseif isempty(Input1) |
---|
40 | [FileName, PathName] = uigetfile('*.mat', 'Input Quadrupole Center file.'); |
---|
41 | if ~isstr(FileName) |
---|
42 | return |
---|
43 | else |
---|
44 | eval(['load ',PathName,FileName]); |
---|
45 | end |
---|
46 | else |
---|
47 | QMS = Input1; |
---|
48 | FileName = []; |
---|
49 | end |
---|
50 | end |
---|
51 | catch |
---|
52 | error('Problem getting input data'); |
---|
53 | end |
---|
54 | |
---|
55 | %QMS = quadplot(QMS); |
---|
56 | |
---|
57 | if nargin < 2 |
---|
58 | FigureHandle = gcf; |
---|
59 | else |
---|
60 | if all(FigureHandle ~= 0) |
---|
61 | figure(FigureHandle); |
---|
62 | end |
---|
63 | end |
---|
64 | |
---|
65 | if all(FigureHandle ~= 0) |
---|
66 | clf reset |
---|
67 | if MontiCarloFlag |
---|
68 | AxesHandles(1) = subplot(3,1,1); |
---|
69 | AxesHandles(2) = subplot(3,1,2); |
---|
70 | AxesHandles(3) = subplot(3,1,3); |
---|
71 | else |
---|
72 | AxesHandles(1) = subplot(2,1,1); |
---|
73 | AxesHandles(2) = subplot(2,1,2); |
---|
74 | end |
---|
75 | |
---|
76 | if QMS.QuadPlane ~= 1 |
---|
77 | set(FigureHandle,'units','normal','position',[Buffer .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]); |
---|
78 | %set(FigureHandle,'units','normal','position',[0+Buffer .09+Buffer .5-2*Buffer .9-2*Buffer-HeightBuffer]); |
---|
79 | else |
---|
80 | set(FigureHandle,'units','normal','position',[.5+.003 .25+Buffer .5-Buffer-.003 .75-1.2*Buffer-HeightBuffer]); |
---|
81 | %set(FigureHandle,'units','normal','position',[.5+Buffer .09+Buffer .5-2*Buffer .9-2*Buffer-HeightBuffer]); |
---|
82 | end |
---|
83 | |
---|
84 | %subplot(3,2,1); |
---|
85 | axes(AxesHandles(1)); |
---|
86 | errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),'b'); |
---|
87 | hold on; |
---|
88 | errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,1), QMS.FitParametersStd(:,1),'.r'); |
---|
89 | errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,1), QMS.FitParametersStd(QMS.GoodIndex,1),'.b'); |
---|
90 | hold off |
---|
91 | title('Linear Fit'); |
---|
92 | ylabel('Y-Intercept'); |
---|
93 | %xlabel('BPM Number'); |
---|
94 | grid |
---|
95 | xaxis([0 length(QMS.FitParameters(:,1))+1]); |
---|
96 | |
---|
97 | |
---|
98 | %subplot(3,1,2); |
---|
99 | axes(AxesHandles(2)); |
---|
100 | errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),'b'); |
---|
101 | hold on |
---|
102 | errorbar(1:length(QMS.FitParameters(:,1)), QMS.FitParameters(:,2), QMS.FitParametersStd(:,2),'.r'); |
---|
103 | errorbar(QMS.GoodIndex, QMS.FitParameters(QMS.GoodIndex,2), QMS.FitParametersStd(QMS.GoodIndex,2),'.b'); |
---|
104 | hold off |
---|
105 | xlabel('BPM Number'); |
---|
106 | ylabel('Slope'); |
---|
107 | %xlabel('BPM Number'); |
---|
108 | grid |
---|
109 | xaxis([0 length(QMS.FitParameters(:,1))+1]); |
---|
110 | end |
---|
111 | |
---|
112 | |
---|
113 | if MontiCarloFlag |
---|
114 | % Monte Carlo the mean and standard deviation of all the good offsets |
---|
115 | N = 20000; |
---|
116 | QMS.OffsetMean = NaN*ones(length(QMS.FitParameters(:,2)),1); |
---|
117 | QMS.OffsetStd = NaN*ones(length(QMS.FitParameters(:,2)),1); |
---|
118 | |
---|
119 | OffsetMean1 = QMS.OffsetMean; |
---|
120 | OffsetStd1 = QMS.OffsetStd; |
---|
121 | NormalRV1 = randn(N,1); |
---|
122 | NormalRV2 = randn(N,1); |
---|
123 | if QMS.QuadraticFit == 1 |
---|
124 | NormalRV3 = randn(N,1); |
---|
125 | end |
---|
126 | for i = QMS.GoodIndex(:)' % 1:length(QMS.Offset) |
---|
127 | rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1; |
---|
128 | rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2; |
---|
129 | if QMS.QuadraticFit == 1 |
---|
130 | rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3; |
---|
131 | x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic); |
---|
132 | else |
---|
133 | x = -rvinter ./ rvslope; |
---|
134 | end |
---|
135 | OffsetMean1(i) = mean(x); |
---|
136 | OffsetStd1(i) = std(x); |
---|
137 | end |
---|
138 | |
---|
139 | OffsetMean2 = QMS.OffsetMean; |
---|
140 | OffsetStd2 = QMS.OffsetStd; |
---|
141 | NormalRV1 = randn(N,1); |
---|
142 | NormalRV2 = randn(N,1); |
---|
143 | if QMS.QuadraticFit == 1 |
---|
144 | NormalRV3 = randn(N,1); |
---|
145 | end |
---|
146 | for i = QMS.GoodIndex(:)' % 1:length(QMS.Offset) |
---|
147 | rvinter = QMS.FitParameters(i,1) + QMS.FitParametersStd(i,1)*NormalRV1; |
---|
148 | rvslope = QMS.FitParameters(i,2) + QMS.FitParametersStd(i,2)*NormalRV2; |
---|
149 | if QMS.QuadraticFit == 1 |
---|
150 | rvquadradic = QMS.FitParameters(i,3) + QMS.FitParametersStd(i,3)*NormalRV3; |
---|
151 | x = (-rvslope + sqrt(rvslope.^2 - 4*rvquadradic.*rvinter)) ./ (2*rvquadradic); |
---|
152 | else |
---|
153 | x = -rvinter ./ rvslope; |
---|
154 | end |
---|
155 | OffsetMean2(i) = mean(x); |
---|
156 | OffsetStd2(i) = std(x); |
---|
157 | end |
---|
158 | |
---|
159 | QMS.OffsetMean(:,1) = [OffsetMean1 + OffsetMean2] / 2; |
---|
160 | QMS.OffsetStd(:,1) = [OffsetStd1 + OffsetStd2] / 2; |
---|
161 | |
---|
162 | |
---|
163 | if all(FigureHandle ~= 0) |
---|
164 | |
---|
165 | % plot([OffsetStd1 OffsetStd2 OffsetStd3],'x') |
---|
166 | |
---|
167 | %subplot(3,1,3); |
---|
168 | axes(AxesHandles(3)); |
---|
169 | %errorbar(1:length(QMS.Offset), QMS.Offset, QMS.OffsetStd,'.r'); |
---|
170 | %hold on |
---|
171 | errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), QMS.OffsetStd(QMS.GoodIndex),'.b'); |
---|
172 | |
---|
173 | |
---|
174 | % errorbar(1:length(QMS.Offset), QMS.Offset, QMS.OffsetStd,'.b'); |
---|
175 | % hold on |
---|
176 | % errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd1,'.r'); |
---|
177 | % errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd2,'.g'); |
---|
178 | % errorbar(1:length(QMS.Offset), QMS.Offset, OffsetStd3,'.k'); |
---|
179 | |
---|
180 | % errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), QMS.OffsetStd(QMS.GoodIndex),'.b'); |
---|
181 | % hold on |
---|
182 | % errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd1(QMS.GoodIndex),'.r'); |
---|
183 | % errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd2(QMS.GoodIndex),'.g'); |
---|
184 | % errorbar(QMS.GoodIndex, QMS.Offset(QMS.GoodIndex), OffsetStd3(QMS.GoodIndex),'.k'); |
---|
185 | |
---|
186 | hold off |
---|
187 | ylabel('X-Intercept = - (Y-Intercept)/Slope'); |
---|
188 | ylabel('Monte Carlo'); |
---|
189 | xlabel('BPM Number'); |
---|
190 | grid |
---|
191 | xaxis([0 length(QMS.Offset)+1]); |
---|
192 | %yaxis([-.025 .025]+QMS.Center); |
---|
193 | end |
---|
194 | |
---|
195 | QMS.OffsetSTDMontiCarlo = sqrt((sum(QMS.OffsetStd(QMS.GoodIndex).^2))/length(QMS.GoodIndex)); |
---|
196 | end |
---|