function out = fit_EffPhase_over_energy_cep(Espec, cep_list, traj_mat, options)
% Calulate effective phase and extrapolation beyond cutoff
% Input: energy spectra axis (Espec [eV]), CEP loop (cep_list [pi]),
% trajectory parameters (traj_mat)
arguments
    Espec
    cep_list
    traj_mat
    options.weight_scale (1,1) double {mustBeNumeric} = 1
end

load_const;
load_color;

Up = traj_mat(1).Up;
N_spec = 100;
W_scale = 1;

out.theta_eff_dir = nan(N_spec, numel(cep_list));
out.Theta_Fit_list_dir = nan(numel(Espec), numel(cep_list));
out.offset_dir = nan(numel(cep_list));
out.Edir = nan(numel(cep_list), 1000);
out.Thetadir = nan(numel(cep_list), 1000);

out.theta_eff_rec = nan(N_spec, numel(cep_list));
out.Theta_Fit_list_rec = nan(numel(Espec), numel(cep_list));
out.offset_rec = nan(numel(cep_list));
out.Erec = nan(numel(cep_list), 1000);
out.Thetarec = nan(numel(cep_list), 1000);

weight_co_scale = 25e-4;

figure(500); subplot(2,1,1); hold off; subplot(2,1,2); hold off

for cep_i = 1:numel(cep_list)
    disp(cep_list(cep_i))

    input = traj_mat(cep_i);
    if isempty(input.Ebirth)
        disp('nothing to see here')
        continue
    end

    roi_x = ~isnan(input.x(1,:)) & input.x(1,:) >= 0;
    max_rate_tot = max(input.rate);

    % direct electrons
    if input.use_dir
        roi_d = input.n_recoll == 0 & roi_x;

        E_curr = input.Ekin(roi_d);
        t_curr = input.t_birth(roi_d);
        rate_curr = input.rate(roi_d);
        S_curr = input.S.T(roi_d) + input.S.V(roi_d) - W_scale*input.S.W(roi_d);

        if numel(t_curr)>2 & max(rate_curr) > 1e-4*max_rate_tot
            weight = rate_curr / max(rate_curr);
            weight = weight.^options.weight_scale;

            cut_off_dir = max(E_curr);

            Espec_hist = linspace(0, 1.1*cut_off_dir, N_spec);
            
            % effective phase
            A_dir = histw_phase(E_curr, rate_curr, S_curr, Espec_hist);
            theta_dir = unwrap(angle(A_dir));
            
            roi_weight = weight > weight_co_scale;
            Etest = max(E_curr(roi_weight));

            roi_smm_dir = Espec_hist < Etest;

            % extrapolation
            if numel(Espec_hist(roi_smm_dir)) < 2
                Theta_fitted = zeros(size(Espec));
            else
                Theta_fitted = interp1(Espec_hist(roi_smm_dir), theta_dir(roi_smm_dir), Espec, "linear", "extrap");
            end

            [~, Eminidx_curr_all] = findpeaks(-abs(E_curr-Etest));
            weight_help = zeros(size(weight)); weight_help(Eminidx_curr_all) = weight(Eminidx_curr_all);
            [~, Eminidx_curr] = max(weight_help);
            [~, Eminidx] = min(abs(Espec-Etest));

            offset = round(S_curr(Eminidx_curr)/(const.hbar*pi) - Theta_fitted(Eminidx)/pi, 0);
            
            subplot(2,1,1)
            scatter(E_curr/input.Up, S_curr/const.hbar/pi, 0.1+20*weight, my_B, 'filled')
            hold on
            plot(Espec/input.Up, Theta_fitted/pi + offset, 'LineWidth', 1.5, 'Color', 'c')
            plot(Espec_hist/input.Up, theta_dir/pi + offset, '--', 'LineWidth', 1, 'Color', 'k')
            xlim([0 2*cut_off_dir/input.Up])
            ylim([-15 15])
            hold off
            drawnow;

            out.theta_eff_dir(:, cep_i) = theta_dir;
            out.Theta_Fit_list_dir(:, cep_i) = Theta_fitted;
            out.offset_dir(cep_i) = offset;
            out.Edir(cep_i,1:numel(E_curr)) = E_curr;
            out.Thetadir(cep_i,1:numel(E_curr)) = S_curr/const.hbar;

        else
            disp('nothing to fit here')
        end
    end

    % recollision electrons
    if input.use_recoll
        roi_r = input.n_recoll == 1 & roi_x;

        E_curr = input.Ekin(roi_r);
        t_curr = input.t_birth(roi_r);
        rate_curr = input.rate(roi_r);
        S_curr = input.S.T(roi_r) + input.S.V(roi_r) - W_scale*input.S.W(roi_r);

        if numel(t_curr)>2 & max(rate_curr) > 1e-4*max_rate_tot
            weight = rate_curr / max(rate_curr);
            weight = weight.^options.weight_scale;

            cut_off_rec = max(E_curr);

            Espec_hist = linspace(0, 1.1*cut_off_rec, N_spec);
            
            % effective phase
            A_rec = histw_phase(E_curr, rate_curr, S_curr, Espec_hist);
            theta_rec = unwrap(angle(A_rec));
            
            roi_weight = weight > weight_co_scale;
            Etest = max(E_curr(roi_weight));

            roi_smm_rec = Espec_hist < Etest;

            % extrapolation
            if numel(Espec_hist(roi_smm_rec)) < 2
                Theta_fitted = zeros(size(Espec));
            else
                Theta_fitted = interp1(Espec_hist(roi_smm_rec), theta_rec(roi_smm_rec), Espec, "linear", "extrap");
            end

            [~, Eminidx_curr_all] = findpeaks(-abs(E_curr-Etest));
            weight_help = zeros(size(weight)); weight_help(Eminidx_curr_all) = weight(Eminidx_curr_all);
            [~, Eminidx_curr] = max(weight_help);
            [~, Eminidx] = min(abs(Espec-Etest));

            offset = round(S_curr(Eminidx_curr)/(const.hbar*pi) - Theta_fitted(Eminidx)/pi, 0);
            
            subplot(2,1,2)
            scatter(E_curr/input.Up, S_curr/const.hbar/pi, 0.1+20*weight, my_R, 'filled')
            hold on
            plot(Espec/input.Up, Theta_fitted/pi + offset, 'LineWidth', 1.5, 'Color', 'm')
            plot(Espec_hist/input.Up, theta_rec/pi+offset, '--', 'LineWidth', 1, 'Color', 'k')
            xlim([0 2*cut_off_rec/input.Up])
            ylim([-5 25])
            drawnow
            hold off

            out.theta_eff_rec(:, cep_i) = theta_rec;
            out.Theta_Fit_list_rec(:, cep_i) = Theta_fitted;
            out.offset_rec(cep_i) = offset;
            out.Erec(cep_i,1:numel(E_curr)) = E_curr;
            out.Thetarec(cep_i,1:numel(E_curr)) = S_curr/const.hbar;

        else
            disp('nothing to fit here')
        end


    end
end

end
