-
Notifications
You must be signed in to change notification settings - Fork 1
/
run_LiTrackOpt_clean.m
153 lines (121 loc) · 4.75 KB
/
run_LiTrackOpt_clean.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
clear all;
fontsize = 14;
show = 1;
% Load spectrum data
load('data_samples_1108.mat');
data_spectrum = spec/sum(spec);
spectrum_axis = spectrum_axis/1000;
% Load wakefield data
global A;
A = load('slac.dat');
% Create Parameter struct
global PARAM;
% Set number of sim steps
ESsteps = 200;
% Initialize parameters
par_lims_1108;
param_1108;
params = zeros(nPar,ESsteps);
pscaled = zeros(nPar,ESsteps);
pCurrent = zeros(nPar,1);
pCurrent(1) = PARAM.INIT.SIGZ0; % Bunch Length
pCurrent(2) = PARAM.INIT.SIGD0; % Initial Energy Spread
pCurrent(3) = PARAM.INIT.NPART; % Number of Particles
pCurrent(4) = PARAM.INIT.ASYM; % Initial Gaussian Asymmetry
pCurrent(5) = PARAM.NRTL.AMPL; % Amplitude of RF Compressor
pCurrent(6) = PARAM.NRTL.PHAS; % RF Compressor Phase
pCurrent(7) = PARAM.NRTL.R56; % RTL compression
pCurrent(8) = PARAM.NRTL.T566; % RTL second order compression
pCurrent(9) = decker; % 2-10 Phase
pCurrent(10) = ramp; % Ramp Phase
pCurrent(11) = PARAM.LI20.ELO; % Low Energy Cutoff
pCurrent(12) = PARAM.LI20.EHI; % High Energy Cutoff
pCurrent(13) = PARAM.LI20.BETA; % Beta Function
pCurrent(14) = PARAM.LI20.R16; % Dispersion
pCurrent(15) = PARAM.LI20.T166; % 2nd Order Dispersion
pCurrent(16) = delta; % Energy offset
% Record initial values
pInit = pCurrent;
% Initialize ES
[w, dt] = init_ES(nPar); % Get frequencies and time step
alpha = 2000; % ES parameter
gain = 2; % ES parameter
cost = zeros(1,ESsteps); % Cost
Part_frac = zeros(1,ESsteps); % Fraction of Particles lost
residual = zeros(1,ESsteps); % Chi2 difference between spectra
if show; figure(1); end;
tic;
for j=1:ESsteps;
disp(j);
% Run LiTrack
OUT = LiTrackOpt('FACETpar');
% Calculate particle fraction
Part_frac(j) = 1 - OUT.I.PART(2)/PARAM.INIT.NESIM;
% Interpolate simulated spectrum
sim_spectrum = interpSimSpec(OUT,spectrum_axis,PARAM.SIMU.BIN,delta,PARAM.LI20.R16);
% Calculate residual
residual(j) = sum(data_spectrum.*(sim_spectrum - data_spectrum).^2);
% Set Cost as the value of the residual + particle fraction
%cost(j) = 20 + log(residual(j)) + 0.001*Part_frac(j);
cost(j) = 20 + log(residual(j));
% ES Calc
pLast = 2*(pCurrent-Cent)./Diff;
pNext = pLast'+dt*cos(w*j*dt+gain*cost(j)).*(alpha*w).^0.5;
pNext(pNext < -1) = -1;
pNext(pNext > 1) = 1;
pCurrent = Diff.*pNext'/2+Cent;
% Update Params
PARAM.INIT.SIGZ0 = pCurrent(1); % Bunch Length
PARAM.INIT.SIGD0 = pCurrent(2); % Initial Energy Spread
PARAM.INIT.NPART = pCurrent(3); % Number of Particles
PARAM.INIT.ASYM = pCurrent(4); % Initial Gaussian Asymmetry
PARAM.NRTL.AMPL = pCurrent(5); % Amplitude of RF Compressor
PARAM.NRTL.PHAS = pCurrent(6); % RF Compressor Phase
PARAM.NRTL.R56 = pCurrent(7); % RTL Compression
PARAM.NRTL.T566 = pCurrent(8); % RTL Second order compression
decker = pCurrent(9); % 2-10 Phase
ramp = pCurrent(10); % Ramp Phase
PARAM.LI20.ELO = pCurrent(11); % Low Energy Cutoff
PARAM.LI20.EHI = pCurrent(12); % High Energy Cutoff
PARAM.LI20.BETA = pCurrent(13); % Beta Function
PARAM.LI20.R16 = pCurrent(14); % Dispersion
PARAM.LI20.T166 = pCurrent(15); % 2nd Order Dispersion
delta = pCurrent(16); % Energy offset
% Record evolving params
params(:,j) = pCurrent;
pscaled(:,j) = pNext;
if show
figure(1);
subplot(2,2,1);
plot(spectrum_axis,data_spectrum,'g',spectrum_axis,sim_spectrum,'b','linewidth',2);
axis([-4 4 0 6e-3]);
xlabel('X (mm)','fontsize',12);
title('Bunch Spectra','fontsize',10);
legend('DATA','SIM');
subplot(2,2,2);
plot(1:ESsteps,residual,'color','r','linewidth',2);
axis([0 ESsteps 0 7e-6]);
xlabel('Step','fontsize',12);
title('Residual','fontsize',10);
subplot(2,2,3);
plot(1:ESsteps,1-Part_frac,'color','g','linewidth',2);
axis([0 ESsteps 0.75 1.1]);
xlabel('Step','fontsize',12);
title('Particle Fraction','fontsize',10);
subplot(2,2,4);
plot(1:ESsteps,cost,'color','b','linewidth',2);
axis([0 ESsteps 0 9]);
xlabel('Step','fontsize',12);
title('Cost','fontsize',10);
end
end
toc;
% Plot Output
if show
figure(2)
plot(spectrum_axis,data_spectrum,'g',spectrum_axis,sim_spectrum,'b');
legend('DATA','ES-SIMULATION');
xlabel('X (mm)','fontsize',14);
text(-3.5,5e-3,['Residual = ' num2str(residual(j),'%0.2e')],'fontsize',14);
stupidest_plotter;
end