function [MaxABRValues] = ABRMaxWaveValues(ee, StimModality)
% anex\ABRMaxWaveValues finds the ABR response with the largest P1 amplitude and
% smallest latency for which there are user W1 annotations
% input: 
%   ee (anex): experiment to investigate
%   StimModality (string): sitmulus modality ('Optical', 'Acoustic') % does
%   not yet work for acoustic!
% output:
% MaxABRValues (struct) with:
%   MaxP1N1Amplitude (double): max P1N1 peak for any presented bera of the
%       sitmModality
%   P1Latency (double) shortest detected latency in any wave
%   AppliedIntensities (double) mW value for the max. wave
%   Hardware (string): used Laser
%   B (berabr) where the max ampl. was found and berabr where the shortest
%   latency was found

MaxABRValues=[]; %return empty if no data exists yet
if status_cache==1
    try
        load_name = strcat(StimModality,'ABRMaxResponse.mat');
        load(fullfile(getProcessedDataDir(ee),load_name)),
        disp('ABR threshold results loaded.')
    catch 
        disp('ABR threshhold results need to be compiled');
    end
elseif status_cache==0
    % find all protocols that contain intensity protocols
    [BoutA,iiTraceOutA] = findModProt(ee, StimModality,'I');
    [BoutB,iiTraceOutB] = findModProt(ee, StimModality,' ');
    % concatinate them
    if isempty(BoutB)
        B_all = BoutA;
        iiTraceOut = iiTraceOutA;
    elseif ~isempty(BoutB)
        B_all = [BoutA,BoutB];
        iiTraceOut = [iiTraceOutA,iiTraceOutB];
    end
    % prepare output strucutre
    MaxABRValues=struct();

    % separate by different Lasers if existing
    stim_types={};
    for ii=1:length(B_all)
        stim_types{ii}=B_all(ii).Stim.stimulusHardware;
    end
    stim_types=unique(stim_types);
    maxAmplValues=[]; % to store max AMptlidue
    P1LatencyValues=[]; % to store max AMptlidue
    Stim_intensities=[];
    all_Bs=[];
    for jj=1:length(stim_types)
        % for all these bera indexes find the minimum intensity and pool them
        IntMaxAmp_temp = nan(numel(B_all),1,1);
        IntLatency_temp = nan(numel(B_all),1,1);
        indexBeraTrace_temp = nan(numel(B_all),1,1);
    
    
        
        for ii = 1 : numel(B_all)
            if ~strcmp(B_all(ii).Stim(1).stimulusHardware,stim_types{jj})
                continue
            end
            [IntMaxAmp_temp(ii),IntLatency_temp(ii),indexBeraTrace_temp(ii)] = findMaxWave(B_all(ii));
        end
        
        % find the lowest intensity
        [maxValue,ixMax] = max(IntMaxAmp_temp); %there might be more than one minimum values, if the same intensity has been used multiple times
        maxAmplValues=[maxAmplValues,maxValue];
        P1LatencyValues=[P1LatencyValues,IntLatency_temp(ixMax)];
        all_Bs=[all_Bs,B_all(ixMax)];
        if ~isnan(maxValue)
            Stim_intensities=[Stim_intensities,B_all(ixMax).C.Ical(indexBeraTrace_temp(ixMax))];
        end

    end

    MaxABRValues.MaxP1N1Amplitude = maxAmplValues;
    MaxABRValues.P1Latency = P1LatencyValues;
    MaxABRValues.AppliedIntensities=Stim_intensities;
    MaxABRValues.Hardware=stim_types;
    MaxABRValues.B = all_Bs;
    save_name = strcat(StimModality,'ABRMaxResponse.mat');
    save(fullfile(getProcessedDataDir(ee),save_name),'MaxABRValues'),
    disp('ABR max Wave results saved')
end
    
end









