source: MML/trunk/machine/SOLEIL/StorageRing/bba_mml/quadcalcoffset.m @ 4

Last change on this file since 4 was 4, checked in by zhangj, 11 years ago

Initial import--MML version from SOLEIL@2013
File size: 11.4 KB
Line 
1function [Offsetx Offsety] = quadcalcoffset(varargin)
2% Compute final BBA offsets, especially for the 40 combined BPMs
3%
4%  INPUTS
5%  Optional
6%  1. 'Write', 'NoWrite' - Create file to be used for application onto Liberas
7%  2. 'Display', 'NoDisplay' - plot computed final offset
8%
9%  See Also quadgetdat, quadplotAll
10
11
12%
13%  NOTES
14%  1. Application
15%     use /home/operateur/GrpDiagnostics/matlab/DserverBPM/Set_BBA_Offsets_planH.m and Set_BBA_Offsets_planV.m
16%     use file tableBBAH.mat and tableBBAV
17%     Set_BBA_Offsets_planH('ADD')
18%     Set_BBA_Offsets_planV('ADD')
19%
20
21%
22% Written by Laurent S. Nadolski
23
24DebugFlag = 0;
25DisplayFlag = 1;
26WriteFlag = 0; % write results for applying to TANGO
27
28for i = length(varargin):-1:1
29    if ~isempty(varargin)
30        if strcmpi(varargin{i}, 'Display')
31            DisplayFlag = 1;
32            varargin(i) = [];
33        elseif strcmpi(varargin{i}, 'NoDisplay')
34            DisplayFlag = 0;
35            varargin(i) = [];
36        elseif strcmpi(varargin{i}, 'Write')
37            WriteFlag = 1;
38            varargin(i) = [];
39        elseif strcmpi(varargin{i}, 'NoWrite')
40            WriteFlag = 0;
41            varargin(i) = [];
42        end
43    end
44end
45
46Offsetx = [];
47Offsety = [];
48
49% List of BPM whose offset is determined by 2 quadrupoles
50doubleBPMListQ2Q3 =[
51    1     3
52    4     7
53    5     3
54    8     7
55    9     3
56    12     7
57    13     3
58    16     7
59    ];
60
61doubleBPMListQ6Q7 = [
62    1     7
63    2     3
64    2     8
65    3     3
66    3     8
67    4     3
68    5     7
69    6     3
70    6     8
71    7     3
72    7     8
73    8     3
74    9     7
75    10     3
76    10     8
77    11     3
78    11     8
79    12     3
80    13     7
81    14     3
82    14     8
83    15     3
84    15     8
85    16     3
86    ];
87
88doubleBPMListQ4Q5=[
89    1     4
90    4     4
91    5     4
92    8     4
93    9     4
94    12     4
95    13     4
96    16     4
97    ];
98
99doubleBPMListQ11Q12=[
100    13    9
101    ];
102
103doubleBPMList = [doubleBPMListQ2Q3 ; doubleBPMListQ4Q5 ; doubleBPMListQ6Q7; doubleBPMListQ11Q12];
104
105% get data from BBA
106[X, Y, BPMxFamily, QUADxFamily, BPMyFamily, QUADyFamily, DateMatX, DateMatY, QUADxDevList, QUADyDevList] = quadgetdata([], 0);
107
108
109DIR0 = pwd;
110cd(getfamilydata('Directory', 'BBAcurrent'));
111
112% Maker saying if offset has been calculated
113% -1 if not enough data
114% -2 if too many files
115%  1 if offset computed
116if ~isempty(X)
117    isDonex=zeros(size(X(:,1)));
118end
119if ~isempty(Y)
120
121    isDoney=zeros(size(Y(:,1)));
122end
123
124if ~isempty(X)
125    % look for double BPM in measdata in H-plane
126    for k=1:size(X(:,1))
127        %isdoubleBPM
128        if isDonex(k) == 0 % offset to be determined
129            if ~isempty(findrowindex(X(k,1:2), doubleBPMList))
130                if DebugFlag
131                    fprintf('BPMx(%d,%d) offset determined by 2 quadrupoles\n', X(k,1), X(k,2));
132                end
133                % look for the 2 set of data on meas data
134                Idx = find(X(:,1) == X(k,1) & X(:,2) == X(k,2));
135                NbOfData = length(Idx);
136                if NbOfData < 2
137                    fprintf('BPMx(%d,%d) only %d data %s(%d,%d): missing data\n', ...
138                        X(k,1), X(k,2), NbOfData, QUADxFamily(k,:), QUADxDevList(k,:));
139                    isDonex(k) = -1;
140                elseif NbOfData > 2
141                    fprintf('BPMx(%d,%d)  %d : too many data\n', X(k,1), X(k,2), NbOfData);
142                    for ik=1:NbOfData,
143                        fprintf('%s(%d,%d)\n', QUADxFamily(Idx(ik),:), QUADxDevList(Idx(ik),:));
144                        str = sprintf(' bba%d%s%dh* ', QUADxDevList(Idx(ik),1), deblank(QUADxFamily(Idx(ik),:)), QUADxDevList(Idx(1),2));
145                        system(['ls -C1 ' str]);
146
147                    end
148                    fprintf('\n')
149                    isDonex(Idx) = -2;
150                else % check data are well from 2 different quadrupoles
151                    if diff(X(Idx,8)) == 0
152                        fprintf('BPMx(%d,%d) not properly determined : same quad %s(%d,%d)\n', X(k,1), X(k,2), QuadxFamily(Idx), X(k,10:11));
153                    else %Compute combined offset
154                        Offsetx = [Offsetx; X(Idx(1),:)];
155                        [tmpOffset tmpSigma]= getCombinedOffset(X(Idx,:), QUADxFamily(Idx,:), QUADxDevList(Idx,:), 1);
156                        Offsetx(end, 3) = tmpOffset;
157                        Offsetx(end, 4) = tmpSigma;
158                        Offsetx(end, 6) = mean(X(Idx,6)); % Mean current
159                        Offsetx(end, 8) = NaN;
160                        Offsetx(end, 9) = NaN;
161                        isDonex(Idx) = 1;
162                    end
163                end
164            else % single BPM-quad offset
165                % Check for too many data
166                Idx = find(X(:,1) == X(k,1) & X(:,2) == X(k,2));
167                NbOfData = length(Idx);
168                if NbOfData > 1
169                    fprintf('BPMx(%d,%d)  too many data (%d)\n', X(k,1), X(k,2), NbOfData);
170                    str = sprintf(' bba%d%s%dh* ', QUADxDevList(Idx(1),1), deblank(QUADxFamily(Idx(1),:)), QUADxDevList(Idx(1),2));
171                    system(['ls -C1 ' str]);
172                    isDonex(Idx) = -2;
173                else
174                    isDonex(k) = 1; % do nothing since computation already done by quadgetdata
175                    Offsetx = [Offsetx; X(k,:)];
176                end
177            end
178        end
179    end
180end
181
182% look for double BPM in measdata in H-plane
183if ~isempty(Y)
184    for k=1:size(Y(:,1))
185        %isdoubleBPM
186        if isDoney(k) == 0 % offset already determined
187            if ~isempty(findrowindex(Y(k,1:2), doubleBPMList))
188                if DebugFlag
189                    fprintf('BPMz(%d,%d) offset determined by 2 quadrupoles\n', Y(k,1), Y(k,2));
190                end
191                % look for the 2 set of data on meas data
192                Idx = find(Y(:,1) == Y(k,1) & Y(:,2) == Y(k,2));
193                NbOfData = length(Idx);
194                if NbOfData < 2
195                    fprintf('BPMz(%d,%d) only %d data %s(%d,%d): missing data\n', ...
196                        Y(k,1), Y(k,2), NbOfData, QUADyFamily(k,:), QUADyDevList(k,:));
197                    isDoney(k) = -1;
198                elseif NbOfData > 2
199                    fprintf('BPMz(%d,%d)  %d: too many data\n', Y(k,1), Y(k,2), NbOfData);
200                    for ik=1:NbOfData,
201                        fprintf('%s(%d,%d)\n', QUADyFamily(Idx(ik),:), QUADyDevList(Idx(ik),:));
202                        str = sprintf(' bba%d%s%dv* ', QUADyDevList(Idx(ik),1), deblank(QUADyFamily(Idx(ik),:)), QUADyDevList(Idx(1),2));
203                        system(['ls -C1 ' str]);
204                    end
205                   
206                    fprintf('\n')
207                    isDoney(Idx) = -2;
208                else % check data are well from 2 different quadrupoles
209                    if diff(Y(Idx,8)) == 0
210                        fprintf('BPMz(%d,%d) not properly determined: same quad %s(%d,%d)\n', Y(k,1), Y(k,2), QUADyFamily(Idx(1),:), QUADyDevList(Idx(1),:));
211                    else %Compute combined offset
212                        Offsety = [Offsety; Y(Idx(1),:)];
213                        [tmpOffset tmpSigma]= getCombinedOffset(Y(Idx,:), QUADyFamily(Idx,:), QUADyDevList(Idx,:), 2);
214                        Offsety(end, 3) = tmpOffset;
215                        Offsety(end, 4) = tmpSigma;
216                        Offsety(end, 6) = mean(Y(Idx,6)); % Mean current
217                        Offsety(end, 8) = NaN;
218                        Offsety(end, 9) = NaN;
219                        isDoney(Idx) = 1;
220                    end
221                end
222            else % single BPM-quad offset
223                % Check for too many data
224                Idx = find(Y(:,1) == Y(k,1) & Y(:,2) == Y(k,2));
225                NbOfData = length(Idx);
226                if NbOfData > 1
227                    fprintf('BPMz(%d,%d)  too many data (%d)\n', Y(k,1), Y(k,2), NbOfData);
228                    str = sprintf(' bba%d%s%dv* ', QUADyDevList(Idx(1),1), deblank(QUADyFamily(Idx(1),:)), QUADyDevList(Idx(1),2));
229                    system(['ls -C1 ' str]);
230                    isDoney(Idx) = -2;
231                else
232                    isDoney(k) = 1; % do nothing since computation already done by quadgetdata
233                    Offsety = [Offsety; Y(k,:)];
234                end
235            end
236        end
237    end
238end
239
240if DisplayFlag
241%     if isempty(Offsetx)
242%         fprintf('Offsetx is empty\n')
243%         return;
244%     end
245%     if isempty(Offsety)
246%         fprintf('Offsety is empty\n')
247%         return;
248%     end
249    % Plot data
250    L = getcircumference;
251
252    figure;
253    subplot(2,1,1);
254    if ~isempty(X)
255        errorbar(Offsetx(:,5), Offsetx(:,3), Offsetx(:,4), '.b');
256    end
257    ylabel('Horizontal [mm]');
258    xaxis([0 L]);
259    title('BPM Offset with BPM determined by 2 quadrupoles');
260
261    subplot(2,1,2);
262    if ~isempty(Y)
263        errorbar(Offsety(:,5), Offsety(:,3), Offsety(:,4), '.b');
264    end
265    xlabel('BPM Position [meters]');
266    ylabel('Vertical [mm]');
267    xaxis([0 L]);
268   
269end
270
271% Write 2 files for applying BBA offsets to TANGO
272
273if WriteFlag
274    if isempty(Offsetx)
275        fprintf('Offsetx is empty\n')
276    else
277        for ik=1:size(Offsetx,1),
278            tableBBAH(ik,1) = dev2tangodev('BPMx', Offsetx(ik,1:2)); tableBBAH{ik,2} = num2str(Offsetx(ik,3), '%+.3f');
279        end
280        save('tableBBAH.mat','tableBBAH')
281    end
282    if isempty(Offsety)
283        fprintf('Offsety is empty\n')
284    else
285        for ik=1:size(Offsety,1),
286            tableBBAV(ik,1) = dev2tangodev('BPMz', Offsety(ik,1:2)); tableBBAV{ik,2} = num2str(Offsety(ik,3), '%+.3f');
287        end
288        save('tableBBAV.mat','tableBBAV')
289    end
290
291    % output result in text format
292    fileName = 'BBAresult.txt';
293    fid = fopen(fileName, 'wt');
294    fprintf(fid, '*************************************************\n');
295    fprintf(fid, '** BBA results:  %s **\n', datestr(now) );
296    fprintf(fid, '*************************************************\n\n');
297    fprintf(fid, '*************** H-plane ******************** \n');
298    if ~isempty(Offsetx)
299        for ih = 1:size(tableBBAH,1)
300            fprintf(fid, '%s # %s  \n', tableBBAH{ih,1:2});
301        end
302    end
303    if ~isempty(Offsety)
304        fprintf(fid, '%s \n', ' ');
305        fprintf(fid, '*************** V-plane ******************** \n');
306        for iv = 1:size(tableBBAV,1)
307            fprintf(fid, '%s # %s  \n', tableBBAV{iv,1:2});
308        end
309    end
310    fclose(fid);
311    system(['nedit ' fileName '&']);
312end
313cd(DIR0);
314
315function [OffsetC SigmaC]= getCombinedOffset(data, Fam, devList, plane)
316%Look for double BPM neigboring two quads
317% plane = 1 H-plane
318% plane = 2 V-plane
319
320DebugFlag = 1;
321
322w = zeros(2,2);
323for k=1:2,
324    w(k,:) = quadgetbetaKL(Fam(k,:), devList(k,:));
325end
326
327OffsetC = (data(1,3)*w(1,plane) + data(2,3)*w(2,plane))/(w(1,plane) + w(2,plane));
328SigmaC = (data(1,4)*w(1,plane) + data(2,4)*w(2,plane))/(w(1,plane) + w(2,plane));
329
330if DebugFlag
331    if plane == 1
332        fprintf('BPMx(%2d,%d) Combined offset: %+f  %s(%2d,%2d) offset: %+f  %s(%2d,%2d) offset: %+f\n', ...
333            data(1,1), data(1,2), OffsetC, Fam(1,:), devList(1,:), data(1,3), Fam(2,:), devList(2,:), data(2,3));
334    else
335        fprintf('BPMz(%2d,%d) Combined offset: %+f  %s(%2d,%2d) offset: %+f  %s(%2d,%2d) offset: %+f\n', ...
336            data(1,1), data(1,2), OffsetC, Fam(1,:), devList(1,:), data(1,3), Fam(2,:), devList(2,:), data(2,3));
337    end
338end
339
340function lambda = getbetaKL(Family, DeviceList)
341
342% Integrated gradient
343KL = getkleff(Family, DeviceList);
344
345lambda = zeros(2,length(KL));
346
347% betax at entrance of magnet
348[bx bz] = modelbeta(Family, DeviceList);
349
350lambda(1,:) = abs(KL).*bx;
351lambda(2,:) = abs(KL).*bz;
Note: See TracBrowser for help on using the repository browser.