%================================================================
%  NLC simple example
%  un-ideal tuned machine, two beamlines
%================================================================
%
clear all;
close all;

% electrons after IP
fname = input('Enter file name of beam1 data: ','s');
beam1      = load(fname);
pbeam=250.;
w0_laser=1.e-3; %1000 micron laser w0
R11=-1.68;
R12=0.;
R21=0.00555;
R22=-0.595; %TRANSPORT from IP to mid-chicane
R33=-2.26;
R34=0.;
R43=-0.099;
R44=-0.443; %TRANSPORT from IP to mid-chicane
i1=[1 0 0 0 0 0];
i2=[0 1 0 0 0 0];
i3=[0 0 1 0 0 0 ];
i4=[0 0 0 1 0 0];
i5=[0 0 0 0 1 0];
i6=[0 0 0 0 0 1];
e_beam1=i1*beam1';
angx_beam1=i2*beam1'*1e6; %angle in urad
angy_beam1=i3*beam1'*1e6; %angle in urad
angx_chicane=R22*angx_beam1;
angy_chicane=R44*angy_beam1;
ang_beam1=sqrt(angx_beam1.^2+angy_beam1.^2);
ang_chicane=sqrt(angx_chicane.^2+angy_chicane.^2);
z_beam1=i4*beam1';
x_beam1=i5*beam1';
y_beam1=i6*beam1';
x_chicane=R11*x_beam1+R12*angx_beam1;
y_chicane=R33*y_beam1+R34*angy_beam1;
DeltaP=(1.-cos((pbeam/0.44065)*ang_beam1*1.e-6));
DeltaP_lumwt=0.25*DeltaP;
DeltaP_chicane=(1.-cos((pbeam/0.44065)*ang_chicane*1.e-6));
r2_beam1=1.e-18*(x_chicane.^2+y_chicane.^2); %units are m^2
laserwt=exp(-2*r2_beam1/w0_laser^2);
DeltaP_chicane_laserwt=(1./sum(laserwt))*laserwt*DeltaP_chicane';
% printout some results on means, stds
mean_e_beam=(mean(e_beam1)-pbeam)/250*1e6 % in ppm
mean_x_beam1=mean(x_beam1)
mean_y_beam1=mean(y_beam1)
mean_z_beam1=mean(z_beam1)
mean_angx_beam1=mean(angx_beam1)
mean_angy_beam1=mean(angy_beam1)
std_e_beam1=std(e_beam1)/pbeam*1e6 % in ppm
std_x_beam1=std(x_beam1)
std_y_beam1=std(y_beam1)
std_z_beam1=std(z_beam1)
std_angx_beam1=std(angx_beam1)
std_angy_beam1=std(angy_beam1)
mean_lumwted_DeltaP=mean(DeltaP_lumwt)
mean_estimated_DeltaP_at_chicane=mean(DeltaP_chicane)
mean_estimated_DeltaP_at_chicane_with_laserwting=mean(DeltaP_chicane_laserwt)

% PLOTS
subplot(2,2,1)
hist(abs(e_beam1),100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Beam1 Energy (GeV)','fontsize',14);
%set(gca,'yscale','log')

subplot(2,2,3)
hist(angx_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,4)
hist(angy_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Y Angle (urad)','fontsize',14);

figure

subplot(2,2,1)
hist(DeltaP_lumwt,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('IP Lum-wted Depol','fontsize',14);

subplot(2,2,2)
plot(sqrt(r2_beam1)*1e3,laserwt,'r.');
grid on
axis([0 5 0 1]);
xlabel('beam radius (mm)','fontsize',14);
ylabel('laser weight','fontsize',14);

subplot(2,2,3)
hist(DeltaP_chicane,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Predicted chicane depol','fontsize',14);

subplot(2,2,4)
hist(DeltaP_chicane_laserwt,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Pred. Chicane depol wted','fontsize',14);

figure
subplot(2,2,1)
hist(x_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron X (nm)','fontsize',14);

subplot(2,2,2)
hist(y_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Y (nm)','fontsize',14);

subplot(2,2,3)
hist(z_beam1,100);
grid on
%axis([0 20 -0.025 0.0]);
xlabel('Electron Z (um)','fontsize',14);

subplot(2,2,4)
plot(z_beam1,e_beam1,'r.')
grid on
xlabel('Electron Z (um)','fontsize',14);
ylabel('Electron Energy (GeV)','fontsize',14);

figure
subplot(2,2,1)
plot(x_beam1,y_beam1,'r.')
grid on
xlabel('Electron X (nm)','fontsize',14);
ylabel('Electron Y (nm)','fontsize',14);

subplot(2,2,2)
plot(angx_beam1,angy_beam1,'r.')
grid on
xlabel('Electron X Angle (urad)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

subplot(2,2,3)
plot(e_beam1,x_beam1,'r.')
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron X (nm)','fontsize',14);

subplot(2,2,4)
plot(e_beam1,y_beam1,'r.')
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron Y (nm)','fontsize',14);

figure
subplot(2,2,1)
plot(e_beam1,angx_beam1,'r.')
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,2)
plot(e_beam1,angy_beam1,'r.')
grid on
xlabel('Electron Energy (GeV)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);

subplot(2,2,3)
plot(x_beam1,angx_beam1,'r.')
grid on
xlabel('Electron X (nm)','fontsize',14);
ylabel('Electron X Angle (urad)','fontsize',14);

subplot(2,2,4)
plot(y_beam1,angy_beam1,'r.')
grid on
xlabel('Electron Y (nm)','fontsize',14);
ylabel('Electron Y Angle (urad)','fontsize',14);







%print -dpdf -f1 beam1_IP_1.pdf
%print -dpdf -f2 beam1_IP_2.pdf
%print -dpdf -f3 beam1_IP_3.pdf
%print -dpdf -f4 beam1_IP_4.pdf
%print -dpdf -f5 beam1_IP_5.pdf

%
%exit