function [x_kl_unique,y_m_square_mean,y_m_square_std,x_kl_fit,y_m_square_fit,chisq_parab,beta_ip_first_meth,sigma_beta_ip_first_meth,beta_ip_second_meth,sigma_beta_ip_second_meth,emittance_ip,sigma_emittance_ip,alpha_ip_first_meth,sigma_alpha_ip_first_meth,alpha_ip_second_meth,sigma_alpha_ip_second_meth]=twiss_analysis_sha(y_m_mean,y_m_std,x_A_unique,emittance_extr_m,sigma_emittance_extr_m,ws_chosen,sig_x_or_y)





%INPUT DATA





%parameters to calculate the twiss parameters



if strcmp(ws_chosen,'WIP')    %at the IP for the optics beta x=0.04 and beta y=0.001

    if strcmp(sig_x_or_y,'X') %for X beam size

      a_par=1.8517915549374;  

    elseif strcmp(sig_x_or_y,'Y')  %for Y beam size

      a_par=1.2364609683258; 

    end

    kl_qd0_design=1.36396800693;   %for alpha estimation

elseif strcmp(ws_chosen,'WPIP Tungsten') || strcmp(ws_chosen,'WPIP Carbon')  %at the PIP for the optics beta x=0.04 and beta y=0.001

    if strcmp(sig_x_or_y,'X')    %for X beam size

      a_par=3.1449945094737;   

    elseif strcmp(sig_x_or_y,'Y') %for Y beam size

      a_par=2.3453424065088; 

    end

    kl_qd0_design=1.112672023613;   %for alpha estimation

end

 







%CONVERSION FOR PARABOLA PARAMETRIZATION



x_kl_unique=x_A_unique*0.0104462;  %conversion A==>KL

y_m_square_mean=y_m_mean.^2; %conversion um==>m^2







%CALCULATION OF THE ERROR SQUARE

y_m_square_std=2*y_m_mean.*y_m_std;







%FITTING



%[A,B,C,rms,chisq] = parabola_fit(x_kl_unique',y_m_square_mean',y_m_square_std');  %without sutracted dispersion

[q_parab,dq_parab,chisq_parab,Cv_parab]=noplot_parab_2(x_kl_unique',y_m_square_mean',y_m_square_std');





%CALCULATION OF TWISS PARAMETERS AT THE IP



beta_ip_first_meth=emittance_extr_m*a_par^2/q_parab(1);      %with first Sha method

beta_ip_second_meth=a_par*sqrt(q_parab(3)/q_parab(1));  %with second Sha method

emittance_ip=sqrt(q_parab(1)*q_parab(3))/a_par;        %with second Sha method

%y_m_square_fit_qd0_design=q_parab(1)*(kl_qd0_design-q_parab(2)).^2+q_parab(3);

%alpha_ip=sqrt((y_m_square_fit_qd0_design/(emittance_extr_m*beta_ip_second_meth))-1);

alpha_ip_first_meth=a_par*(q_parab(2)-kl_qd0_design)/beta_ip_first_meth;

alpha_ip_second_meth=a_par*(q_parab(2)-kl_qd0_design)/beta_ip_second_meth;









%CALCULATION OF THE ERRORS OF THE TWISS PARAMETERS



sigma_beta_ip_first_meth=a_par^2*sqrt((1/q_parab(1))^2*sigma_emittance_extr_m^2+(emittance_extr_m/q_parab(1)^2)^2*dq_parab(1)^2);

sigma_emittance_ip=1/a_par*sqrt((q_parab(3)/(2*sqrt(q_parab(1)*q_parab(3))))^2*dq_parab(1)^2 + (q_parab(1)/(2*sqrt(q_parab(1)*q_parab(3))))^2*dq_parab(3)^2);

%sigma_emittance_ip_cov=1/a_par*sqrt((q_parab(3)/(2*sqrt(q_parab(1)*q_parab(3))))^2*dq_parab(1)^2 + (q_parab(1)/(2*sqrt(q_parab(1)*q_parab(3))))^2*dq_parab(3)^2+2*q_parab(3)/(2*sqrt(q_parab(1)*q_parab(3)))*q_parab(1)/(2*sqrt(q_parab(1)*q_parab(3)))*Cv_parab(1,3));



sigma_beta_ip_second_meth=a_par*sqrt((q_parab(3)/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1)^2))^2*dq_parab(1)^2 + (1/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1)))^2*dq_parab(3)^2);

%sigma_beta_ip_second_meth_cov=a_par*sqrt((q_parab(3)/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1)^2))^2*dq_parab(1)^2 + (1/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1)))^2*dq_parab(3)^2-2*q_parab(3)/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1)^2)*1/(2*sqrt(q_parab(3)/q_parab(1))*q_parab(1))*Cv_parab(1,3));



%sigma_y_m_square_fit_qd0_design=sqrt((kl_qd0_design-q_parab(2))^4*dq_parab(1)^2 + (2*q_parab(2)*q_parab(1)-2*q_parab(1)*kl_qd0_design)^2*dq_parab(2)^2 + dq_parab(3)^2);

%sigma_alpha_ip= sqrt((1/(2*beta_ip_second_meth*emittance_extr_m*sqrt((y_m_square_fit_qd0_design-beta_ip_second_meth*emittance_extr_m)/(beta_ip_second_meth*emittance_extr_m))))^2*sigma_y_m_square_fit_qd0_design^2+(y_m_square_fit_qd0_design/(2*beta_ip_second_meth*emittance_extr_m^2*sqrt((y_m_square_fit_qd0_design-beta_ip_second_meth*emittance_extr_m)/(beta_ip_second_meth*emittance_extr_m))))^2*sigma_emittance_extr_m^2+(y_m_square_fit_qd0_design/(2*beta_ip_second_meth^2*emittance_extr_m*sqrt((y_m_square_fit_qd0_design-beta_ip_second_meth*emittance_extr_m)/(beta_ip_second_meth*emittance_extr_m))))^2*sigma_beta_ip_second_meth^2);

sigma_alpha_ip_first_meth= sqrt((a_par*(x_kl_fit_waist-kl_qd0_design)/beta_ip_first_meth^2)^2*sigma_beta_ip_first_meth^2);

sigma_alpha_ip_second_meth= sqrt((a_par*(x_kl_fit_waist-kl_qd0_design)/beta_ip_second_meth^2)^2*sigma_beta_ip_second_meth^2);







end