forked from ThorstenHellert/SC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SCgetBPMreading.m
276 lines (236 loc) · 7.36 KB
/
SCgetBPMreading.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
function [B,T] = SCgetBPMreading(SC,varargin)
% SCgetBPMreading
% ===============
%
% NAME
% ----
% SCgetBPMreading - Calculates BPM readings based on the current injection scheme.
%
% SYNOPSIS
% --------
% `[B, T] = SCgetBPMreading(SC [, options])`
%
%
% DESCRIPTION
% -----------
% This function calculates the particle trajectories based on the current
% injection setup as defined in `SC.INJ` and calculates the corresponding BPM
% readings. The injection setup is a structure with the fields:
%
% `nParticles`::
% Number of particles for tracking.
% `nTurns`::
% Number of turns for tracking.
% `nShots`::
% Number of injections used for averaging the BPM readings.
% `Z0ideal`::
% [6 x 1] array defining the ideal injection vector.
% `Z0`::
% [6 x 1] array defining the current injection vector.
% `beamSize`::
% [6 x 6] array defining the beam sigma matrix.
% `randomInjectionZ`::
% [6 x 1] array defining the shot-to-shot injection jitter.
% `trackMode`::
% String defining the tracking mode. If set to orbit mode ('ORB'), the AT
% function `findorbit6` is used to calculate the trajectories. Otherwise,
% bunches are generated and tracking is performed. In both cases the
% corresponding BPM readings are calculated. If the tracking mode is
% 'pORB', the pseudo-orbit is calculated by averaging the turn-by-turn
% BPM readings.
%
% If the global variable `plotFunctionFlag` is 1, the tracking results are plotted.
%
%
% INPUTS
% ------
%
% `SC`::
% SC base structure
%
% OPTIONS
% -------
% The following options can be given as name/value-pairs:
%
% `'BPMords'` (`SC.ORD.BPM`):: List of BPM ordinates at which the reading should be returned
%
% RETURN VALUES
% -------------
% `B`::
% BPM readings
% `T`::
% Particle trajectories
%
% GLOBALS
% -------
% `plotFunctionFlag`:: If true, each BPM reading is plotted
%
% EXAMPLES
% --------
% Switch on plotting, track 100 particles for 2 turns and store the BPM readings in `B`.
% ------------------------------------------------------------------
% global plotFunctionFlag
% plotFunctionFlag = 1;
% SC.INJ.trackMode = 'TBT';
% SC.INJ.nParticles = 100;
% SC.INJ.nTurns = 2;
% B = SCgetBPMreading(SC);
% ------------------------------------------------------------------
%
% SEE ALSO
% --------
% *SCgenBunches*, *SCregisterBPMs*, *SCplotBPMreading*
global plotFunctionFlag
% Parse optional arguments
p = inputParser;
addOptional(p,'BPMords',[]);
parse(p,varargin{:});
% Switch orbit or tracking mode
if strcmp(SC.INJ.trackMode,'ORB')
nTurns = 1;
nParticles = 1;
else
nTurns = SC.INJ.nTurns;
nParticles = SC.INJ.nParticles;
end
% Prelocate BPM readings
B1 = nan(2,nTurns*length(SC.ORD.BPM),SC.INJ.nShots);
% Check if plotting is desired
if plotFunctionFlag
T1 = nan(6,nTurns*nParticles*length(SC.RING),SC.INJ.nShots);
refOrds = 1:length(SC.RING);
else
refOrds = SC.ORD.BPM;
end
% Loop over injections
for nShot=1:SC.INJ.nShots
% Switch orbit or tracking mode
if strcmp(SC.INJ.trackMode,'ORB')
% Calculate closed orbit
T = findorbit6(SC.RING, refOrds);%,SC.INJ.Z0);
else
% Generate initial particle distribution
Zin = SCgenBunches(SC);
% Calculate trajectories
T = atpass(SC.RING, Zin, 1, nTurns, refOrds);
end
% Account for particle lost (AT sets only x to nan)
T(:,isnan(T(1,:))) = nan;
if plotFunctionFlag
% Store for plotting
T1(:,:,nShot) = T;
end
% Get BPM reading
B1(:,:,nShot) = calcBPMreading(SC,T,'atAllElements',plotFunctionFlag);
end
% Calculate average BPM reading
B = mean(B1,3,'omitnan');
% Plot trajectories
if plotFunctionFlag
SCplotBPMreading(SC,B,T1);
end
% Get pseudo-orbit BPM reading
if strcmp(SC.INJ.trackMode,'pORB')
Bpseudo = nan(2,length(SC.ORD.BPM));
for nBPM=1:length(SC.ORD.BPM)
Bpseudo(:,nBPM) = mean(B(:,nBPM:length(SC.ORD.BPM):end),2,'omitnan');
end
B = Bpseudo;
end
% Check if user specified list of BPMs has been defined
if ~isempty(p.Results.BPMords)
% Select BPMs for return array
ind = find(ismember(SC.ORD.BPM,p.Results.BPMords));
% Check if all ordinates have been recognized
if length(ind)~=length(p.Results.BPMords)
warning('Not all specified ordinates are registered BPMs.')
end
% Get multiturn indices
if strcmp(SC.INJ.trackMode,'TBT')
ind = [0:(nTurns-1)]*length(SC.ORD.BPM)+ind';
end
% Define final array of BPM readings
B = B(:,ind(:));
end
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary function
function B = calcBPMreading(SC,T,varargin)
% Calculates the BPM readings based on a multi-{particle,turn,shot} trajectories. The beam is
% considered lost if the relative amount of lost particles exceeds the value specified in
% `SC.INJ.beamLostAt` and the corresponding BPM reading is `NaN`.
%
% INPUTS
% ------
% `SC`:: SC base structure.
% `T`:: Particle trajectories, e.g. as given by `atpass`.
%
% OPTIONS
% -------
% 'atAllElements' (0)::
% Flags if input trajectories are given only at BPM ordinates or at all lattice elements.
%
% RETURN VALUES
% -------------
% `B`::
% [2 x (length(BPMords)*nTurns)] array of x and y BPM readings [m]
% Parse optional arguments
p = inputParser;
addOptional(p,'atAllElements',0);
parse(p,varargin{:});
% Switch orbit or tracking mode
if strcmp(SC.INJ.trackMode,'ORB')
nTurns = 1;
BPMnoise = cell2mat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'NoiseCO'))';
nParticles = 1;
else
nTurns = SC.INJ.nTurns;
BPMnoise = cell2mat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'Noise'))';
nParticles = SC.INJ.nParticles;
end
% Define BPM errors along trajectory (multiple turns)
BPMoffset = repmat(cell2mat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'Offset'))' + cell2mat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'SupportOffset'))',1,nTurns);
BPMcalError = repmat(cell2mat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'CalError'))',1,nTurns);
BPMroll = repmat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'Roll')',1,nTurns) + repmat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'SupportRoll')',1,nTurns);
BPMnoise = repmat( BPMnoise ,1,nTurns) .* SCrandnc(2,2,nTurns*length(SC.ORD.BPM));
BPMsumError = repmat(atgetfieldvalues(SC.RING(SC.ORD.BPM),'SumError')',1,nTurns);
% Check if trajectories were given not only at BPMs
if p.Results.atAllElements
% Indices of BPMs in trajectories
nE = reshape((0:nTurns-1)*length(SC.RING)+SC.ORD.BPM',1,[]);
else
nE = 1:(length(SC.ORD.BPM)*nTurns);
end
% Check for multiparticle tracking
if nParticles > 1
% Reshape trajectories in tensor
M = SCparticlesIn3D(T,nParticles);
% Read x and y positions
Tx = squeeze(M(1,nE,:));
Ty = squeeze(M(3,nE,:));
% Calculate center of charge
Bx1 = mean(Tx,2,'omitnan')';
By1 = mean(Ty,2,'omitnan')';
% Check for detactable BPM signal
beamLost = find( sum( isnan( Tx ),2 )' .* (1 + BPMsumError.*SCrandnc(2,size(BPMsumError))) > ( nParticles * SC.INJ.beamLostAt ),1);
% Reflect effective non detectable signal
Bx1(beamLost:end) = nan;
By1(beamLost:end) = nan;
else
% Read x and y positions
Bx1 = T(1,nE);
By1 = T(3,nE);
end
% Add roll error
Bx = cos(BPMroll) .* Bx1 - sin(BPMroll) .* By1;
By = sin(BPMroll) .* Bx1 + cos(BPMroll) .* By1;
% Add calibration and offset
Bx = (Bx - BPMoffset(1,:)) .* (1+BPMcalError(1,:));
By = (By - BPMoffset(2,:)) .* (1+BPMcalError(2,:));
% Add noise
Bx = Bx + BPMnoise(1,:);
By = By + BPMnoise(2,:);
% Store BPM reading in output variable
B(1,:) = Bx;
B(2,:) = By;
end