function patchfile_analysis(trial,ddata,patchlog_name,hy,stim_artifact_threshold,stim_artifact_jitter,IPI_jitter,raw_window,baseline_time_window,PSC_time_window,Q_time_window,...
        Q_btwn_thresholds,PSCbegin_threshold,PSCend_threshold,asynchronous_window,mcc_output_voltage_scaling_factor,samples_per_ms,permutations_of_IPIs,number_of_IPIs,possible_repetitions_of_permutation,...
        APs_in_the_train,max_length_raw,min_length_raw,metrics_of_interest,measured_ITI,real_time_mode,amillion,atrillion,samples_per_second,number_of_trials,plot_controls,save_that_sht,...
        asynch_window_fractionofsecond,minimum_possible_decay_time,maximum_possible_decay_time,user_vetting,do_careful_latency_measuring,do_detecting,do_fitting,do_riseplateaumeasurement,do_montecarlojitteredevent,...
        threshforfitting,fitendtime,successful_sims_required,maximum_possible_simulation_loops,do_smoothing,do_multipeak_detecting,display_traces,save_after_each_trial,...
        do_asynch_carefuldetection,do_I_Na_checking,do_PSCclass_measuring,do_controls_measuring,do_raw_measuring,do_access_input_measuring,do_fitting_to_fast_and_slow_templates,fastandslow_criteria_window)
    
%plot raw recording alone
if ~real_time_mode
    figure(1)
    plot(ddata)
    ylim([-2 1])
    title(strcat('Raw recording__',patchlog_name,'__',num2str(trial)))
    set(1,'Position',[16 263 1575 501])
end

%these will become handy below when saving
asky='-ascii';dub='-double';tbs='-tabs';shortpatch_name=strrep(patchlog_name,'.hdf',''); 
shortpatch_name=strrep(shortpatch_name,'H:\work_modulebayHD\','');
common_pathname=strcat('H:\extracted_values_modulebayHD\',shortpatch_name,'__');

%extracting the experimental info stored at the beginning of each recording
OVSF_presyn=ddata(1,1);   %the same information could be had from either column
sizeddata=size(ddata);
post_oct_19_2003=1;
post_feb_11_2004=1;
if post_oct_19_2003         
    OVSF_postsyn=ddata(2,1);
    pulse_duration_presyn=ddata(3,1);
    pulse_duration_postsyn=ddata(4,1);
    sampling_rate=ddata(5,1);
    stim_philosophy=ddata(6,1);             %can set at 6 for forced reads of hy=[2 1;50 50;100 1] type stimuli as hy=[25 50 100] type stimuli
    stim_philosophy=6   %MUST remove
    output_channel_to_presyn=ddata(7,1);
    ddata=ddata(8:end,:);
    if post_feb_11_2004
        temperature_during_train_mean=ddata(1,1)*10; %multiply to account for conversion factor of the heater/thermometer (.1V/degree C)
        temperature_during_train_std=ddata(2,1)*10;
        ddata=ddata(3:end,:);
    else
        temperature_during_train_mean=NaN;
        temperature_during_train_std=NaN;
    end
else
    pulse_duration_presyn=1;
    OVSF_postsyn=NaN;
    pulse_duration_postsyn=NaN;
    sampling_rate=NaN;
    output_channel_to_presyn=NaN;
    temperature_during_train_mean=NaN;
    temperature_during_train_std=NaN;
    ddata=ddata(2:end,:);
end
shrunksizeddata=size(ddata);
%ddata(:,1)=ddata(:,2);%MUST remove -- for measuring autaptic currents

%converting the recorded signal, voltage from the board, into pA, by the conversion factor used by multiclamp commander
ddata=ddata*mcc_output_voltage_scaling_factor;
% if trial<2, ddata(:,1)=ddata(:,1)*10; end %MUST CHANGE BACK -- this is bc 2 different gain settings on mcc were in use (added in this code 9/27/2006)

stim_artifact_jitter=stim_artifact_jitter*(samples_per_ms*pulse_duration_presyn);%used in finding the stims from the recorded trains to compensate for temporal variability
%anything within this value following the peak of a pulse artifact will be considered part of that same pulse

leak(1)=mean(ddata(1:1000,1));ddata(:,1)=ddata(:,1)-leak(1); %subtracting leak, so thresholds etc. will work properly
leak(2)=mean(ddata(1:1000,2));ddata(:,2)=ddata(:,2)-leak(2); 

raitpd=100;   %duration of the Raccess/input test pulse used during the expt        %MUST CHANGE BACK to look at older expts with shorter test pulses                                        
for j=1:2   %finding stim times then measuring input and access resistance for each channel
   stim_peaks=find_peaks_within_threshold_crossing_regions(stim_artifact_threshold,stim_artifact_jitter,ddata(:,j))   %find stim peaks
   %if any(find(trial==[43]))&&j==1,stim_peaks=[30000 stim_peaks];end%MUST remove%manual set event times
   stim_indices(1:length(stim_peaks),j)=stim_peaks'  
    
   %find the access resistance and input resistance
   if do_access_input_measuring
       [r_acc_min_val(j), r_acc_min_ind(j), access_resistance(j), input_resistance(j)]=find_access_and_input_resistances_from_test_pulse(ddata(:,j),stim_peaks(1),raitpd,samples_per_ms);
   end
end

if do_careful_latency_measuring %for this peak detection, unlike the general one used above, the goal is to find the end of the peak;  only chan 1 is considered, and the threshold should be high enough to miss the R_acc test pulse
    %the actual measurment of the template-match criterion that will be used for the measure of latency is done below;  here is for detection of the ends of the peaks (onset of I_Na)
    presyn_INa_onsets=find_peaks_specifically_ends_of_peaks(.9*max(ddata(:,1)),50,ddata(:,1))
%     if any(find(trial==[1:100])),
%         figure,plot(ddata(:,1)),hold on, plot(presyn_INa_onsets,.9*max(ddata(:,1)),'go'),xlim([presyn_INa_onsets(1)-40 presyn_INa_onsets(1)+40])%plotting
%         figure,plot(ddata(:,1)),hold on, plot(presyn_INa_onsets,.9*max(ddata(:,1)),'go'),xlim([presyn_INa_onsets(2)-40 presyn_INa_onsets(2)+40])
%         figure,plot(ddata(:,1)),hold on, plot(presyn_INa_onsets,.9*max(ddata(:,1)),'go'),xlim([presyn_INa_onsets(3)-40 presyn_INa_onsets(3)+40])
%     end
    %if find(trial==[1]),presyn_INa_onsets=[50974 54975 55475];end %MUST REMOVE -- manual override for detection
    if length(presyn_INa_onsets)~=3, error('wrong number of ends of presyn I_Na onsets detected'),end
end

%checking to make sure the expected number of pulses is present
if stim_philosophy==1
	if ~isempty(hy)
        if length(stim_indices(:,1))~=APs_in_the_train+1  %add one for the hyperpolarizing pulse
            error('see here1')
        end
	else
        if length(stim_indices(:,1))~=2     %one for the real stim artifact, one for the hyperpolarizing test pulse
            error ('see here3')
        end
	end
elseif stim_philosophy==2
    if test_pulse_present & (length(stim_indices(:,1))~=11) %NOT BACK COMPATIBLE
        error('see here4')
    elseif ~test_pulse_present & (length(stim_indices(:,1))~=10)
        error('see here5')
    end
elseif stim_philosophy==3
    %any number of stimuli goes  
elseif stim_philosophy==6   
    if length(stim_indices(:,1))~=4 error('see here6'), end        %counting the test pulse, 4 stimuli
end

%removing the test pulse stim from the list and finding the interpulse intervals
test_pulse_end_time=stim_indices(1,1);
stim_indices=stim_indices(2:end,:) %prevents the hyperpolarizing seal test pulse from being treated as a stimulus pulse (NOT BACK COMPATIBLE)
stim_indices(find(~stim_indices))=NaN %convert the zeros to NaNs
%find the IPI values in the output train to presynaptic
IPI_val_vec=diff(stim_indices(:,1));
IPI_val_vec=IPI_val_vec/samples_per_ms  %converting into ms

%determining the permutation of IPIs presented on this trial
placefinder=[];
if stim_philosophy~=6
	if ~isempty(hy)
        for jia=1:length(hy(1,:))
            for yi=1:hy(1,jia)
                if abs(IPI_val_vec(jia)-(hy(2,jia)+((yi-1)*hy(3,jia))))<IPI_jitter
                    placefinder(jia)=yi;
                    break
                end
            end
		end
        azaz=placefinder
	else
        azaz=[0]
    end
else    %if stimphil is 6 
    azaz=hy(find(IPI_val_vec(1)<(hy+IPI_jitter)&IPI_val_vec(1)>(hy-IPI_jitter)));
    azaz=[azaz, IPI_val_vec(2)]
end

ddata(:,2)=ddata(:,2)-mean(ddata(stim_indices(1,1)+[-200*samples_per_ms:-5*samples_per_ms],2));   %re-correcting for the leak immediately preceding the first stim

%checking for inward current in the presyn trace during the voltage step evoking each AP, to confirm I_Na is present
if do_I_Na_checking
    persistent presyn_I_Na_array
    mama=size(presyn_I_Na_array)
    [presyn_I_Na_present_thistrial]=presyn_I_Na_checker(stim_indices, ddata(:,1),pulse_duration_presyn*samples_per_ms,trial);
    presyn_I_Na_array=[presyn_I_Na_array, presyn_I_Na_present_thistrial];
    
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'presyn_I_Na_array');save(pathname,'presyn_I_Na_array',asky,dub,tbs)
    end
end

%measuring the template-match criterion from the onset of the presyn I_Na (or at least the best estimate thereof) for use in latency measurement
if do_careful_latency_measuring
    persistent criteria_for_latency
    criteria_for_latency_this_trial=[];
    %fastemplate=monoexprisedecay_template_generator([14.16371820548567  11.31185487577630  19.35826590564763 0 1],1:50);fastemplate=[0 fastemplate];
    fastemplate=lsqfit_template_generator([3 8 52 225 0 0],[1:50]);fastemplate=[0 fastemplate]; %this template is much more accurate, at least for some cases (eg 20060413) than the above
    for captain_caveman=1:length(presyn_INa_onsets)
            %removing the artifact from the end of the presyn stim
        deadspace=(pulse_duration_presyn*samples_per_ms)+stim_indices(captain_caveman,1)+[-3:3];
        a=ddata(deadspace(1),2);b=ddata(deadspace(2),2);if a==b,stepper=repmat(a,7,1);else stepper=[a:(b-a)/6:b]';end
        %tempdata=ddata(:,2);  tempdata(1:(deadspace(1)-1))=nan;tempdata((deadspace(end)+1):end)=nan;  before=tempdata(presyn_INa_onsets(captain_caveman)+[0:199]);
        ddata(deadspace,2)=stepper;
        %after=ddata(presyn_INa_onsets(captain_caveman)+[0:199],2);
        %figure, plot(before,'r'),hold on, plot(after,'b')
        postsyn_trace_after_this_AP=ddata(presyn_INa_onsets(captain_caveman)+[0:199],2)';   %using 1 to 199 permits detection of events beginning up to 14.9 ms after the onset of presyn I_Na, assuming 5 ms template
        criterion_this_AP=template_based_mini_detector(postsyn_trace_after_this_AP,fastemplate)';
        criteria_for_latency_this_trial=[criteria_for_latency_this_trial;criterion_this_AP];
        figure,plot(postsyn_trace_after_this_AP,'k'),hold on,plot(criterion_this_AP,'r'), plot([1 199],[2 2],'b'),drawnow
    end
    criteria_for_latency=[criteria_for_latency criteria_for_latency_this_trial];
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'criteria_for_latency');save(pathname,'criteria_for_latency',asky,dub,tbs)
    end
end

%event detection -- ideally this should be done with a template having parameters (rise, decay, etc.) best for the synapse currently being studied
if do_detecting
    persistent criteriaz
	criteriaz_this_trial=[];
	templatte=lsqfit_template_generator([3 8 52 225 0 0],[1:100]);%templatte=[0 templatte]; %it's unclear whether better to include zero at the beginning
	for ayo=1:length(stim_indices(:,1))
        %trace=ddata(stim_indices(ayo,1)+[(pulse_duration_presyn+1)*samples_per_ms:(pulse_duration_presyn+1+PSC_time_window(2))*samples_per_ms],2)';
        trace=ddata(stim_indices(ayo,1)+[20:(PSC_time_window(end)+length(templatte))],2)';
        critt=template_based_mini_detector(trace,templatte)';
        %figure(20+trial),subplot(1,3,ayo),plot(trace,'k'),hold on,plot(critt,'r'), ylim([-5 5])
        criteriaz_this_trial=[criteriaz_this_trial;critt];
	end
	criteriaz=[criteriaz criteriaz_this_trial];
	if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'criterion_vals');save(pathname,'criteriaz',asky,dub,tbs)
    end
end

%curve fitting
if do_fitting
	persistent fittingz
	fittingz_this_trial=[];
	for dibaa=1:length(stim_indices(:,1))   %need to go into iterative_curve_fit to change things like latency etc.
        critlength=length(critt);   %this assumes do_detecting is 1
        currentcrit=criteriaz_this_trial(((dibaa-1)*critlength)+[1:critlength]); 
          if  find(currentcrit > threshforfitting) 
              exceeder=find(currentcrit>threshforfitting);
              exceeder=exceeder(1)-10;
              tracce=ddata(stim_indices(dibaa,1)+exceeder:stim_indices(dibaa,1)+exceeder+fitendtime,2);%figure(100+dibaa),hold on, plot(tracce,'g')
              %trace=ddata(stim_indices(dibaa,1)+[(pulse_duration_presyn+1)*samples_per_ms:21*samples_per_ms],2);       %use this for fitting without having detected
              templatte=iterative_curve_fit(tracce)';
              %figure(50+dibaa),plot(tracce,'k'),hold on, plot(templatte,'r')
          else
              templatte=repmat(nan,fitendtime+12,1);
          end
        fittingz_this_trial=[fittingz_this_trial;templatte];
	end
	fittingz=[fittingz fittingz_this_trial];
	if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'fit_vals');save(pathname,'fittingz',asky,dub,tbs)
    end
end

if do_riseplateaumeasurement
    persistent plateauz
    plateauz_this_trial=[];
    for doowop=1:length(stim_indices(:,1))
        critlength=length(critt);   %this assumes do_detecting is 1
        currentcrit=criteriaz_this_trial(((doowop-1)*critlength)+[1:critlength]);
        if find(currentcrit > threshforfitting)
            exceeder=find(currentcrit>threshforfitting);
            exceeder=exceeder(1)-10;
            raw_segment=ddata(stim_indices(doowop,1)+exceeder:stim_indices(doowop,1)+exceeder+fitendtime,2); %figure(100+doowop),hold on, plot(raw_segment,'r')
            start_point=ceil(fittingz_this_trial((fitendtime+7)+((doowop-1)*(fitendtime+12))))
            peak_point=ceil(fittingz_this_trial((fitendtime+11)+((doowop-1)*(fitendtime+12))))+start_point-1
            [plateau_onset, plateau_length]=find_riseplateau(start_point, peak_point, raw_segment)
            %figure(60+doowop),plot(raw_segment,'k'),hold on, plot(plateau_onset+[0:plateau_length-1],repmat(raw_segment(plateau_onset),1,plateau_length),'go')
        else
            plateau_onset=nan
            plateau_length=nan
        end
        plateauz_this_trial=[plateauz_this_trial;plateau_onset;plateau_length];
    end
    plateauz=[plateauz,plateauz_this_trial];
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'plateau_vals');save(pathname,'plateauz',asky,dub,tbs)
    end
end

if do_montecarlojitteredevent
    persistent montejitterz
    montejitterz_this_trial=[];
    for bebop=1:length(stim_indices(:,1))
        critlength=length(critt);    %this assumes do_detecting is 1
        currentcrit=criteriaz_this_trial(((bebop-1)*critlength)+[1:critlength]);
        fitted_minimal_psc=fittingz_this_trial((fitendtime+8)+((bebop-1)*(fitendtime+12)))
        if any(currentcrit > threshforfitting) && (fitted_minimal_psc < -1)
           stored_jitterz=montecarlo_jittered_events(fitted_minimal_psc,plateauz_this_trial(bebop*2),successful_sims_required,maximum_possible_simulation_loops)'
        else 
           stored_jitterz=repmat(nan,successful_sims_required+1,1)
        end
     montejitterz_this_trial=[montejitterz_this_trial;stored_jitterz];
    end
    montejitterz=[montejitterz montejitterz_this_trial]
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'montejittervals_vals');save(pathname,'montejitterz',asky,dub,tbs)
    end
end

%multipeak detection
if do_multipeak_detecting
    persistent multipeakz
    multipeakz_this_trial=[stim_indices(1,1)];
    numstims=length(stim_indices(:,1));
    for illuvitar=2:numstims
        trace_start=stim_indices(illuvitar,1)+samples_per_ms*(pulse_duration_presyn+1);
        if illuvitar<numstims, trace_end=stim_indices(illuvitar+1,1)-(5*samples_per_ms);else trace_end=stim_indices(illuvitar,1)+(300*samples_per_ms);end
        datatrace=ddata(trace_start:trace_end,2);
        holder=find_multiple_peaks(datatrace);
        holder(1:5:end)=holder(1:5:end)+stim_indices(illuvitar,1)+(pulse_duration_presyn+1)*samples_per_ms;    %making all times 'absolute'
        holder(3:5:end)=holder(3:5:end)+stim_indices(illuvitar,1)+(pulse_duration_presyn+1)*samples_per_ms;
        multipeakz_this_trial=[multipeakz_this_trial;holder];
    end
    figure(2),clf,plot(ddata(:,2)),hold on
    avav=multipeakz_this_trial(2:end)
    plot(avav(1:5:end),avav(2:5:end),'ro','markersize',15,'linewidth',3)
    plot(avav(3:5:end),avav(5:5:end),'gs','markersize',15,'linewidth',3)
    error('s')
% adddd=size(multipeakz)
% bdddd=size(multipeakz_this_trial)   
    if size(multipeakz_this_trial,1)>size(multipeakz,1)      %correcting size, bc the number of multipeaks can't be predicted in advance
        multipeakz=cat(1,multipeakz,repmat(nan,(size(multipeakz_this_trial,1)-size(multipeakz,1)),size(multipeakz,2)));
    elseif size(multipeakz,1)>size(multipeakz_this_trial,1)
        multipeakz_this_trial=cat(1,multipeakz_this_trial,repmat(nan,(size(multipeakz,1)-size(multipeakz_this_trial,1)),1));
    end
    multipeakz=[multipeakz multipeakz_this_trial];
   
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'multipeak_vals');save(pathname,'multipeakz',asky,dub,tbs)
	end
end
   
%the approach is to run the template-match for event detection over the specified window, if criterion exceeds thresh, do a curve fit, keep the vals
%do_asynch_carefuldetection=[do it or not, begin window, end window, criterion, samples preceding peak to include in fitting]
if do_asynch_carefuldetection(1)
    persistent asynch_criteriaz
    templatte=lsqfit_template_generator([3 8 52 225 0 0],[1:60]);%templatte=[0 templatte]; %it's unclear whether better to include zero at the beginning
	prestimulation_asynch_trace=ddata(test_pulse_end_time+(150*samples_per_ms):test_pulse_end_time+(1500*samples_per_ms),2)';
    prestimulation_criteriaz_this_trial=template_based_mini_detector(prestimulation_asynch_trace,templatte)';
    %figure(33)
    %plot(prestimulation_asynch_trace,'k'), hold on, plot(prestimulation_criteriaz_this_trial,'r'),error('t')
    prestimulation_asynch_trace=[inf prestimulation_asynch_trace(1:length(prestimulation_criteriaz_this_trial)) inf];
    prestimulation_criteriaz_this_trial=[inf; prestimulation_criteriaz_this_trial; inf];    
    
    asynch_trace=ddata(stim_indices(end,1)+[do_asynch_carefuldetection(2)*samples_per_ms:do_asynch_carefuldetection(3)*samples_per_ms],2)';
    %figure,plot(asynch_trace),ggg=size(asynch_trace)
    asynch_criteriaz_this_trial=template_based_mini_detector(asynch_trace,templatte)';
    %figure,plot(asynch_criteriaz_this_trial),mmm=size(asynch_criteriaz_this_trial)
    asynch_criteriaz=[asynch_criteriaz [[prestimulation_asynch_trace asynch_trace(1:length(asynch_criteriaz_this_trial))]' [prestimulation_criteriaz_this_trial; asynch_criteriaz_this_trial]]];
    %figure(55),plot(asynch_trace,'k'),hold on
    %plot(asynch_criteriaz,'r')
    
%     persistent asynch_fittingz
%     asynch_event_peaks=find_peaks_within_threshold_crossing_regions(do_asynch_carefuldetection(4),5*samples_per_ms,asynch_criteriaz_this_trial)
%     asynch_fittingz_this_trial=[trial;do_asynch_carefuldetection(2)];
%     if asynch_event_peaks
%         longer_asynch_trace=ddata(stim_indices(end,1)+[do_asynch_carefuldetection(2)*samples_per_ms:(do_asynch_carefuldetection(3)+20)*samples_per_ms],2)';
%         for i=1:length(asynch_event_peaks)
%             fittable_asynch=longer_asynch_trace((asynch_event_peaks(i)-do_asynch_carefuldetection(5)):(asynch_event_peaks(i)+200-1-do_asynch_carefuldetection(5)));
%             %the actual time the event happened is latency found by the fit + 3d AP time + beginning of detection window - do_asynch_carefuldetection(5) VERIFY -- could be off by sample
%             asynch_fit=iterative_curve_fit(fittable_asynch');
%             %plot([asynch_event_peaks(i)-do_asynch_carefuldetection(5):asynch_event_peaks(i)-do_asynch_carefuldetection(5)+99],asynch_fit(1:100),'b')
%             %vvvzv=size(asynch_fit)
%             asynch_fittingz_this_trial=[asynch_fittingz_this_trial;asynch_fit'];
% 		end                
%     end
%     if length(asynch_fittingz_this_trial)>size(asynch_fittingz,1)      %correcting size, bc the number of multipeaks can't be predicted in advance
%         asynch_fittingz=cat(1,asynch_fittingz,repmat(nan,(length(asynch_fittingz_this_trial)-size(asynch_fittingz,1)),size(asynch_fittingz,2)));
%     elseif size(asynch_fittingz,1)>length(asynch_fittingz_this_trial)
%         asynch_fittingz_this_trial=cat(1,asynch_fittingz_this_trial,repmat(nan,(size(asynch_fittingz,1)-length(asynch_fittingz_this_trial)),1));
%     end     
%     asynch_fittingz=[asynch_fittingz asynch_fittingz_this_trial]; 
        
	if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'asynch_criterion_vals');save(pathname,'asynch_criteriaz',asky,dub,tbs)
        %pathname=strcat(common_pathname,'asynch_fit_vals');save(pathname,'asynch_fittingz',asky,dub,tbs)
	end
end

if do_fitting_to_fast_and_slow_templates
    if ~exist('fastemplate')
        persistent fastandslow_criteriaz
        %persistent fastemplate,persistent slotemplate
        fastemplate=monoexprisedecay_template_generator([14.16371820548567  11.31185487577630  19.35826590564763 0 1],1:50);    %amplitude, rise, decay, offset, latency;  the first 3 are average of a (small, fast) preinduction and a similar postinduction event;  a latency of 1 is used to prefix the waveform with a zero
        slotemplate=monoexprisedecay_template_generator([10.431819085631   12.507638858663   159.893180756647 0 1],1:50);           %from the slow event after the second AP, trial 5, patch20070319T125639
    end
    datatrace=ddata(stim_indices(1)+fastandslow_criteria_window,2);
    fastandslow_criteriaz_this_trial=fitting_to_fast_and_slow_templates(datatrace',fastemplate,slotemplate);
    %figure,plot(datatrace),'k'),hold on,plot(fastandslow_this_trial(:,1),'b'),plot(fastandslow_this_trial(:,2),'r'),error('t')
    fastandslow_criteriaz=[fastandslow_criteriaz datatrace(1:size(fastandslow_criteriaz_this_trial,1)) fastandslow_criteriaz_this_trial];
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
        pathname=strcat(common_pathname,'fastandslow_criterion_vals');save(pathname,'fastandslow_criteriaz',asky,dub,tbs)
	end    
end


%smoothing -- this is a boxcar smooth, also wrote gaussianwindowsmooth but that is slow (there may be a built-in function for that that's better)
if do_smoothing(1)
	smoothing_binwidth=do_smoothing(2);
	smoothness=movingwindowsmooth(ddata(:,2),smoothing_binwidth);
	ddata(floor((smoothing_binwidth-1)/2):floor((smoothing_binwidth-1)/2)+length(smoothness)-1,2)=smoothness';  %want the smoothed waveform to have the best possible temporal alignment with the original
end

if do_PSCclass_measuring
	%this loop is going to extract all 6 parameters following each AP in each relevant one of the 4 possible synapses
	stsd_parameters={'PSC','charge','latency','risetime','decaytime','asynchronous'};
	stimulated_recorded={{'pre-' 'post-'} {'pre' 'post'}};
	persistent ooverall_values
	for stimulated=1:1%length(stimulated_recorded{1})
         for recorded=2:2%1:length(stimulated_recorded{2})
             %getting the baseline immediately before the first AP stimulation
             if ~isnan(stim_indices(1,stimulated))
                 orig_baseline(recorded)=mean(ddata(stim_indices(1,stimulated)+baseline_time_window,recorded))
                 baselineSD=std(ddata(stim_indices(1,stimulated)+baseline_time_window,recorded));
                 starterthresh=orig_baseline(recorded)-PSCbegin_threshold*baselineSD  %for area and such
                 enderthresh=orig_baseline(recorded)-PSCend_threshold*baselineSD
                 for jb=1:length(find(~isnan(stim_indices(:,stimulated))))      %fill in any additional spaces w/ NaN
                     index_now=stim_indices(jb,stimulated);%below, index_now is sometimes used and stim_indices(jb,stimulated) is also used;  there is some rhyme/reason to this
                     PSC_minz=index_now+PSC_time_window(1)-1+find(ddata(index_now+[PSC_time_window],recorded)==min(ddata(index_now+[PSC_time_window],recorded)))
                     PSC_minindex(jb,stimulated,recorded)=PSC_minz(1)
                     if do_smoothing(1)        %no need to average peak PSC to over nearby points to elim noise when smoothing has already been done
                         PSC(jb,stimulated,recorded)=ddata(PSC_minindex(jb,stimulated,recorded),recorded)-orig_baseline(recorded)
                     else                     
                         PSC(jb,stimulated,recorded)=mean(ddata([PSC_minindex(jb,stimulated,recorded)-(round(.3*samples_per_ms)):PSC_minindex(jb,stimulated,recorded)...
                                 +(round(1*samples_per_ms))],recorded))-orig_baseline(recorded)
                     end
                     if stimulated==recorded,index_now=index_now+2*samples_per_ms;end%to give breathing space when measuring response on the same chan stimulated
                     PSC_startoff=find(ddata(index_now+[PSC_time_window],recorded)<starterthresh); %for finding area (charge) etc.
                     if ~isempty(PSC_startoff)
                        PSC_start_time(jb,stimulated,recorded)=PSC_startoff(1)+index_now+PSC_time_window(1)-1;
                        PSC_endup=find(ddata(PSC_minindex(jb,stimulated,recorded)+minimum_possible_decay_time:PSC_minindex(jb,stimulated,recorded)+maximum_possible_decay_time,recorded)>enderthresh);
                        if ~isempty(PSC_endup)
                             PSC_end_time(jb,stimulated,recorded)=PSC_endup(1)+PSC_minindex(jb,stimulated,recorded)+minimum_possible_decay_time-1;
                             if PSC_end_time(jb,stimulated,recorded)<PSC_minindex(jb,stimulated,recorded),error('PSC ends too early, 194'),end %jaybee=jb,stimm=stimulated,recc=recorded,aa=PSC_end_time,bb=PSC_minindex,end %
                             if Q_btwn_thresholds
                                charge(jb,stimulated,recorded)=sum(ddata(PSC_start_time(jb,stimulated,recorded):PSC_end_time(jb,stimulated,recorded),recorded))
                                psc_duration(jb,stimulated,recorded)=PSC_end_time(jb,stimulated,recorded)-PSC_start_time(jb,stimulated,recorded)+1;%used in plotting
                                charge_window_fractionofsecond(jb,stimulated,recorded)=psc_duration(jb,stimulated,recorded)/samples_per_second                               
                             else
                                charge(jb,stimulated,recorded)=sum(ddata(index_now+[Q_time_window],recorded));
                                charge_window_fractionofsecond(jb,stimulated,recorded)=length(Q_time_window)/samples_per_second
                             end
                            charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)/atrillion; %convert from pA to A
                            charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)*charge_window_fractionofsecond(jb,stimulated,recorded) %correcting for the window length to get coulombs
                            %if charge(jb,stimulated,recorded)>0, charge(jb,stimulated,recorded)=0,end
                            latency(jb,stimulated,recorded)=(PSC_start_time(jb,stimulated,recorded)-stim_indices(jb,stimulated))/samples_per_ms;
                            risetime(jb,stimulated,recorded)=(PSC_minindex(jb,stimulated,recorded)-PSC_start_time(jb,stimulated,recorded))/samples_per_ms;
                            decaytime(jb,stimulated,recorded)=(PSC_end_time(jb,stimulated,recorded)-PSC_minindex(jb,stimulated,recorded))/samples_per_ms;
                        else    %for times when PSC starts but doesn't end
                            disp ('PSC does not end')
                                charge(jb,stimulated,recorded)=sum(ddata(index_now+[Q_time_window],recorded));
                                charge_window_fractionofsecond(jb,stimulated,recorded)=length(Q_time_window)/samples_per_second
                                charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)/atrillion; %convert from pA to A
                                charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)*charge_window_fractionofsecond(jb,stimulated,recorded) %correcting for the window length to get coulombs
                            latency(jb,stimulated,recorded)=NaN;
                            risetime(jb,stimulated,recorded)=NaN;
                            decaytime(jb,stimulated,recorded)=NaN;
                        end
                    else    %when no PSC starts
                        disp ('NO PSC DETECTED')
                        %charge(jb,stimulated,recorded)=NaN;
                                charge(jb,stimulated,recorded)=sum(ddata(index_now+[Q_time_window],recorded));
                                charge_window_fractionofsecond(jb,stimulated,recorded)=length(Q_time_window)/samples_per_second
                                charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)/atrillion; %convert from pA to A
                                charge(jb,stimulated,recorded)=charge(jb,stimulated,recorded)*charge_window_fractionofsecond(jb,stimulated,recorded) %correcting for the window length to get coulombs
                        latency(jb,stimulated,recorded)=NaN;
                        risetime(jb,stimulated,recorded)=NaN;
                        decaytime(jb,stimulated,recorded)=NaN;
                    end
                        if jb==4,    %only reading asynchronous after the fourth pulse (switched from 3d to 4th pulse feb 11, 2007)
                                asynchronous(jb,stimulated,recorded)=sum(ddata(stim_indices(jb,stimulated)+[asynchronous_window],recorded))
                                asynchronous(jb,stimulated,recorded)=asynchronous(jb,stimulated,recorded)/atrillion;%picoAmps to Amps   
                                asynchronous(jb,stimulated,recorded)=asynchronous(jb,stimulated,recorded)*asynch_window_fractionofsecond%to coulombs
                        else
                            asynchronous(jb,stimulated,recorded)=nan;asynchronousforplotting=nan;
                        end
                        %if asynchronous(jb,stimulated,recorded)>0,asynchronous(jb,stimulated,recorded)=0,end
                                
                    ooverall_values(number_of_IPIs+jb,trial,1,stimulated,recorded)=-1*PSC(jb,stimulated,recorded);     %now (oct 10, 2004) using positive PSC, charge and asynch
                    ooverall_values(number_of_IPIs+jb,trial,2,stimulated,recorded)=-1*charge(jb,stimulated,recorded);
                    ooverall_values(number_of_IPIs+jb,trial,3,stimulated,recorded)=latency(jb,stimulated,recorded);
                    ooverall_values(number_of_IPIs+jb,trial,4,stimulated,recorded)=risetime(jb,stimulated,recorded);
                    ooverall_values(number_of_IPIs+jb,trial,5,stimulated,recorded)=decaytime(jb,stimulated,recorded);
                    ooverall_values(number_of_IPIs+jb,trial,6,stimulated,recorded)=asynchronous(jb,stimulated,recorded);
                     
                    ooverall_values(1:number_of_IPIs,trial,1:6,stimulated,recorded)=cat(3,azaz',azaz',azaz',azaz',azaz',azaz');
                end
             end
             %displaying the values (eventually, st they can be approved, but for now, just for observation)
             if display_traces
				figure(2)
				clf
				title(strcat('Trial #',num2str(trial)))
				hold on
				xlabel('Sample'),ylabel('Membrane current (pA)')
				plot(ddata(:,1),'b','linewidth',3)
				plot(ddata(:,2),'k','linewidth',3)
				plot(r_acc_min_ind(1),r_acc_min_val(1),'rd','linewidth',4)
				plot(r_acc_min_ind(2),r_acc_min_val(2),'gd','linewidth',4)
                plot(stim_indices(:,1),ddata(stim_indices(:,1),1),'ro','LineWidth',4)
				plot(PSC_minindex(:,1,2),ddata(PSC_minindex(:,1,2),2),'r+','markersize',10,'LineWidth',4)
				rectangle('position',[stim_indices(1,1)+baseline_time_window(1) -baselineSD baseline_time_window(end)-baseline_time_window(1) 2*baselineSD],'LineWidth',2,'EdgeColor',[.6 0 .6]),hold on
				latprepost=ooverall_values(number_of_IPIs+1:end,trial,3,1,2);risprepost=ooverall_values(number_of_IPIs+1:end,trial,4,1,2);decayprepost=ooverall_values(number_of_IPIs+1:end,trial,5,1,2);
				latprepost=latprepost*samples_per_ms;risprepost=risprepost*samples_per_ms;decayprepost=decayprepost*samples_per_ms;
				if ~find(isnan(latprepost)) %if any (not all, right?) of the values are non-NaN
                    plot(stim_indices(:,1)+latprepost,ddata((stim_indices(:,1)+latprepost),2),'rs','LineWidth',2)
                    plot(stim_indices(:,1)+latprepost+risprepost,ddata(stim_indices(:,1)+latprepost+risprepost,2),'rp','LineWidth',2)
                    plot(stim_indices(:,1)+latprepost+risprepost+decayprepost,ddata(stim_indices(:,1)+latprepost+risprepost+decayprepost,2),'rh','LineWidth',2)
				end
				for presynstim=1:length(stim_indices(:,1))
                    if ooverall_values(number_of_IPIs+presynstim,trial,2,1,2)>0  %if a charge has been detected/stored
                        scaled_charge_prepost=(ooverall_values(number_of_IPIs+presynstim,trial,2,1,2))/charge_window_fractionofsecond(presynstim,1,2)
                        if Q_btwn_thresholds, heighty=scaled_charge_prepost/psc_duration(presynstim,1,2);
                            rectangle('position',[PSC_start_time(presynstim,1,2) -heighty psc_duration(presynstim,1,2) heighty],'LineWidth',2, 'EdgeColor','r')
                        else heighty=scaled_charge_prepost/(Q_time_window(end)-Q_time_window(1)); 
                            rectangle('position',[stim_indices(presynstim,1)+Q_time_window(1) -heighty (Q_time_window(end)-Q_time_window(1)) heighty],'LineWidth',2, 'EdgeColor','r')
                        end
                    end
                    if ooverall_values(number_of_IPIs+presynstim,trial,6,1,2)>0 %if asynchronous release has been detected/stored
		%                 aa=number_of_IPIs
		%                 bb=presynstim
		%                 cc=trial
		%                 dd=ooverall_values(number_of_IPIs+presynstim,trial,6,1,2)
                        heighti=((ooverall_values(number_of_IPIs+presynstim,trial,6,1,2))/asynch_window_fractionofsecond)/(asynchronous_window(end)-asynchronous_window(1));%area divided by length gives height
                        heighti=500*heighti;%scaling for display salience
                        rectangle('position',[stim_indices(presynstim,1)+asynchronous_window(1) -heighti asynchronous_window(end)-asynchronous_window(1) heighti],'LineWidth',2, 'EdgeColor','r')
                    end    
				end
                    
				second_col=stim_indices(:,2)
				for postsynstim=1:length(second_col(~isnan(second_col)))
                   a=PSC_minindex
                   plot(stim_indices(postsynstim,2),ddata(stim_indices(postsynstim,2),2),'go','LineWidth',2)
                   if ~isnan(PSC_minindex(postsynstim,2,1))
                       plot(PSC_minindex(postsynstim,2,1),ddata(PSC_minindex(postsynstim,2,1),1),'g+','LineWidth',2)
                       latpostpre=ooverall_values(number_of_IPIs+1:end,trial,3,2,1),rispostpre=ooverall_values(number_of_IPIs+1:end,trial,4,2,1),decaypostpre=ooverall_values(number_of_IPIs+1:end,trial,5,2,1)
                       latpostpre=latpostpre*samples_per_ms;rispostpre=rispostpre*samples_per_ms;decaypostpre=decaypostpre*samples_per_ms;
                       plot(stim_indices(postsynstim,2)+latpostpre,ddata((stim_indices(postsynstim,2)+latpostpre),1),'gs','LineWidth',2)
                       plot(stim_indices(postsynstim,2)+latpostpre+rispostpre,ddata(stim_indices(postsynstim,2)+latpostpre+rispostpre,1),'gp','LineWidth',2)
                       plot(stim_indices(postsynstim,2)+latpostpre+rispostpre+decaypostpre,ddata(stim_indices(postsynstim,2)+latpostpre+rispostpre+decaypostpre,1),'gh','LineWidth',2)
                       if ooverall_values(number_of_IPIs+postsynstim,trial,2,2,1)>0
                           scaled_charge_postpre=(ooverall_values(number_of_IPIs+postsynstim,trial,2,2,1))/charge_window_fractionofsecond(postsynstim,2,1)
                           if Q_btwn_thresholds,heighty_post=scaled_charge_postpre/psc_duration(postsynstim,2,1); 
                               rectangle('position',[PSC_start_time(postsynstim,2,1) -heighty_post psc_duration(postsynstim,2,1) heighty_post],'linewidth',2,'edgecolor','g')
                           else heighty_post=scaled_charge_postpre/(Q_time_window(end)-Q_time_window(1));
                               rectangle('position',[stim_indices(postsynstim,2)+Q_time_window(1) -heighty_post (Q_time_window(end)-Q_time_window(1)) heighty_post],'linewidth',2,'edgecolor','g')
                           end
                       end
                       if ooverall_values(number_of_IPIs+postsynstim,trial,6,2,1)>0
                           heighti_post=((ooverall_values(number_of_IPIs+postsynstim,trial,6,2,1))/asynch_window_fractionofsecond)/(asynchronous_window(end)-asynchronous_window(1)); %area divided by length gives height
                           heighti=500*heighti;
                           rectangle('position',[stim_indices(postsynstim,2)+asynchronous_window(1) -heighti_post asynchronous_window(end)-asynchronous_window(1) heighti_post],'LineWidth',2, 'EdgeColor','g')
                       end   
                   end
				end
					
			        set(2,'position',[47 694 560 420]),two_obj=flipud(findobj(2));   %flipping so that things in the next 3 plots will plot in the right sequence
                    ylimmm=[-60 9];
			        figure(8),clf,title(strcat('Trial #',num2str(trial))),xlabel('Sample'),ylabel('Membrane current (pA)')
			            h8=subplot(1,1,1);copyobj(two_obj(1:end-2),h8),set(8,'position',[981 694 560 420]),xlim([stim_indices(1,1)-500 stim_indices(2,1)+500]),ylim(ylimmm)
			        figure(9),clf,title(strcat('Trial #',num2str(trial))),xlabel('Sample'),ylabel('Membrane current (pA)')
			            h9=subplot(1,1,1);copyobj(two_obj(1:end-2),h9),set(9,'position',[47 65 481 420]),xlim([stim_indices(1,1)-500 stim_indices(1,1)+500]),ylim(ylimmm)
			        figure(10),clf,title(strcat('Trial #',num2str(trial))),xlabel('Sample'),ylabel('Membrane current (pA)')
			            h10=subplot(1,1,1);copyobj(two_obj(1:end-2),h10),set(10,'position',[576 65 481 420]),xlim([stim_indices(2,1)-500 stim_indices(2,1)+500]),ylim(ylimmm)
			            if size(stim_indices,1)>2
			            figure(11),clf,title(strcat('Trial #',num2str(trial))),xlabel('Sample'),ylabel('Membrane current (pA)')
			                h11=subplot(1,1,1);copyobj(two_obj(1:end-2),h11),set(11,'position',[1099 65 481 420]),xlim([stim_indices(3,1)-500 stim_indices(3,1)+500]),ylim(ylimmm)
			            end
              
                if user_vetting
					user_accepted=input('Do these measurements pass muster? (1 for yes, 0 for no)')
					if ~user_accepted, 
                        ooverall_values(number_of_IPIs+1,trial,1,stimulated,recorded)=input('What is PSC_o?  (Use positive value.)')
                        ooverall_values(number_of_IPIs+2,trial,1,stimulated,recorded)=input('What is PSC_a?  (Use positive value.)')
					end %how to make this work optimally, i.e., how can the user input this information by pointing and clicking on the plot? -- probably need a GUI
                else
                    drawnow
                    pause(.1)
                end
            end 
             if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht)
                 for parameter=1:length(stsd_parameters)        %save the 6 parameters (for the current one of the four possible synapses) separately
                     filename=strcat(stimulated_recorded{1}(stimulated),stimulated_recorded{2}(recorded),stsd_parameters{parameter});
                     pathname=strcat(common_pathname,filename);
                     single_param=ooverall_values(1:end,1:end,parameter,stimulated,recorded);
                     var='single_param';
					 save(pathname{1},var,asky,dub,tbs)
				end
            end
         end
     end
 end
 
%bump counting module
%  bump_indices_and_values=bump_counter(data_with_bumps,samples_per_ms,smoothing_binwidth,network_event_detection_mode,consecutive_ups,consecutive_downs,...
%     asynchronous_event_detection_mode,event_detection_threshold,event_jitter,orig_baseline_postsynaptic,bump_skip_portion,plot_bumps,delete_artifactual_event,patchlog_name,...
%     trial,pulse_num)
% persistent bump_indices_holder, persistent bump_values_holder
% if trial==1, bump_indices_holder=[];bump_values_holder=[];end
% for sweep_segment=1:4
%     if sweep_segment==1,data_with_bumps=ddata(1:(test_pulse_end_time-(102*samples_per_ms)),recorded);adder=0;
%     elseif sweep_segment==2,data_with_bumps=ddata(test_pulse_end_time+samples_per_ms:stim_indices(1)-samples_per_ms,recorded);adder=test_pulse_end_time+samples_per_ms;
%     elseif sweep_segment==3,data_with_bumps=ddata(stim_indices(1)+(3*samples_per_ms):stim_indices(2)-samples_per_ms,recorded);adder=stim_indices(1)+(3*samples_per_ms);
%     elseif sweep_segment==4,data_with_bumps==ddata(stim_indices(2)+(3*samples_per_ms):end,recorded);adder=stim_indices(2)+(3*samples_per_ms);
%     end
%     bump_indices_and_values=bump_counter(data_with_bumps,10,20,1,9,4,0,-3,2,0,1,1,1,'joe',trial,1)
%     error('e')
%     indices=indices+adder;
%     bump_indices_holder(sweep_segment,trial)=1;
%     bump_values_holder(sweep_segment,trial)=1;
% end
                     
if do_controls_measuring
    %a list of the variables, stored row-by-row, in the control figure (columns represent trials;  everything is on one sheet)
	trials=trial;
	ITI=measured_ITI;
	temperature=temperature_during_train_mean; 
	intensity_pre=OVSF_presyn; 
	pulsedur_pre=pulse_duration_presyn; 
	r_access_chanone=access_resistance(1);  %the method of filling these 3, and their equivalents for the other channel, may change
	r_input_chanone=input_resistance(1); 
	leak_chanone=leak(1); 
	tot_charge_chanone=sum(ddata(1:end,1)); 
	intensity_post=OVSF_postsyn;
	pulsedur_post=pulse_duration_postsyn;
	r_access_chantwo=access_resistance(2); 
	r_input_chantwo=input_resistance(2); 
	leak_chantwo=leak(2); 
	tot_charge_chantwo=sum(ddata(1:end),2);
	psc_o_prepre=ooverall_values(number_of_IPIs+1,trial,1,1,1); 
	q_o_prepre=ooverall_values(number_of_IPIs+1,trial,2,1,1);
	lat_prepre=ooverall_values(number_of_IPIs+1,trial,3,1,1);
	ris_prepre=ooverall_values(number_of_IPIs+1,trial,4,1,1);
	decay_prepre=ooverall_values(number_of_IPIs+1,trial,5,1,1);
	asynch_prepre=ooverall_values(number_of_IPIs+1,trial,6,1,1);
	psc_o_prepost=ooverall_values(number_of_IPIs+1,trial,1,1,2); 
	q_o_prepost=ooverall_values(number_of_IPIs+1,trial,2,1,2);
	lat_prepost=ooverall_values(number_of_IPIs+1,trial,3,1,2);
	ris_prepost=ooverall_values(number_of_IPIs+1,trial,4,1,2);
	decay_prepost=ooverall_values(number_of_IPIs+1,trial,5,1,2);
	asynch_prepost=ooverall_values(number_of_IPIs+1,trial,6,1,2);
	if size(ooverall_values,4)>1 %if this analysis includes trials where the postsyn cell was stimulated to fire an AP
		psc_o_postpre=ooverall_values(number_of_IPIs+1,trial,1,2,1);
        q_o_postpre=ooverall_values(number_of_IPIs+1,trial,2,2,1); 
        lat_postpre=ooverall_values(number_of_IPIs+1,trial,3,2,1);
		ris_postpre=ooverall_values(number_of_IPIs+1,trial,4,2,1);
		decay_postpre=ooverall_values(number_of_IPIs+1,trial,5,2,1);
		asynch_postpre=ooverall_values(number_of_IPIs+1,trial,6,2,1);
		psc_o_postpost=ooverall_values(number_of_IPIs+1,trial,1,2,2); 
        q_o_postpost=ooverall_values(number_of_IPIs+1,trial,2,2,2);
        lat_postpost=ooverall_values(number_of_IPIs+1,trial,3,2,2);
		ris_postpost=ooverall_values(number_of_IPIs+1,trial,4,2,2);
		decay_postpost=ooverall_values(number_of_IPIs+1,trial,5,2,2);
		asynch_postpost=ooverall_values(number_of_IPIs+1,trial,6,2,2);
	else
        psc_o_postpre=NaN;q_o_postpre=NaN;lat_postpre=NaN;ris_postpre=NaN;decay_postpre=NaN;asynch_postpre=NaN;psc_o_postpost=NaN;q_o_postpost=NaN;
        lat_postpost=NaN;ris_postpost=NaN;decay_postpost=NaN;asynch_postpost=NaN;
	end
	
	control_variables=[trials ITI temperature intensity_pre pulsedur_pre r_access_chanone r_input_chanone leak_chanone tot_charge_chanone intensity_post ...
            pulsedur_post r_access_chantwo r_input_chantwo leak_chantwo tot_charge_chantwo psc_o_prepre q_o_prepre lat_prepre ris_prepre decay_prepre ...
            asynch_prepre psc_o_prepost q_o_prepost lat_prepost ris_prepost decay_prepost asynch_prepost psc_o_postpre q_o_postpre lat_postpre ris_postpre...
            decay_postpre asynch_postpre psc_o_postpost q_o_postpost lat_postpost ris_postpost decay_postpost asynch_postpost];
	
	control_variable_list={'trials' 'ITI' 'temperature' 'intensity_pre' 'pulsedur_pre' 'r_access_chanone' 'r_input_chanone' 'leak_chanone' 'tot_charge_chanone' 'intensity_post' ...
            'pulsedur_post' 'r_access_chantwo' 'r_input_chantwo' 'leak_chantwo' 'tot_charge_chantwo' 'psc_o_prepre' 'q_o_prepre' 'lat_prepre' 'ris_prepre' 'decay_prepre' 'asynch_prepre'...
            'psc_o_prepost' 'q_o_prepost' 'lat_prepost' 'ris_prepost' 'decay_prepost' 'asynch_prepost' 'psc_o_postpre' 'q_o_postpre' 'lat_postpre' 'ris_postpre' 'decay_postpre'...
            'asynch_postpre' 'psc_o_postpost' 'q_o_postpost' 'lat_postpost' 'ris_postpost' 'decay_postpost' 'asynch_postpost'};
	persistent controlll
	if trial==1, controlll=[];end
	controlll(:,trial)=control_variables';
	
    if (trial==number_of_trials & save_that_sht) | (save_after_each_trial & save_that_sht),filename=strcat(common_pathname,'controls'); var='controlll';save(filename,var,asky,dub,tbs),end
	if trial==number_of_trials&plot_controls,figure(5);set(5,'PaperOrientation','landscape','PaperPosition',[0.5 0.5 10 7.5]);
		hold on
		x_factors={'trials'
		%         'ITI_history'
        %          'temperature'
		%         'intensity_pre' 
		%         'pulsedur_pre'
                  'r_access_chanone'
                  'r_input_chanone'
                  'leak_chanone'
		%         'intensity_post'
		%         'pulsedur_post'
                  'r_access_chantwo'
                  'r_input_chantwo'
                  'leak_chantwo'
                   };
		x_cols=size(x_factors,1);
		y_factors={
        %          'psc_o_prepre'
        %          'q_o_prepre'
        %          'lat_prepre'
        %          'ris_prepre'
        %          'decay_prepre'
        %          'asynch_prepre'
                  'psc_o_prepost'
                  'q_o_prepost'
                  'lat_prepost'
                  'ris_prepost'
                  'decay_prepost'
                  'asynch_prepost'
        %          'psc_o_postpre'
        %          'q_o_postpre'
        %          'lat_postpre'
        %          'ris_postpre'
        %          'decay_postpre'
        %          'asynch_postpre'
        %          'psc_o_postpost'
        %          'q_o_postpost'
        %          'lat_postpost'
        %          'ris_postpost'
        %          'decay_postpost'
        %          'asynch_postpost'
        %          'tot_charge_chanone'
        %          'leak_chanone'
        %          'r_access_chanone'
        %          'r_input_chanone'
        %          'tot_charge_chantwo'
        %          'leak_chantwo'
        %          'r_access_chantwo'
        %          'r_input_chantwo'
        %          'temperature'
                    };
		y_rows=size(y_factors,1);
		for j=1:y_rows  %from convert_stored_matrix just add one more loop to represent epoch, and change the color of the symbol for each epoch
           for k=1:x_cols
               subplot(y_rows,x_cols,((j-1)*x_cols)+k)
               jj=y_factors{j};
               kk=x_factors{k};
               jj=strcmp(y_factors(j),control_variable_list);
               kk=strcmp(x_factors(k),control_variable_list);
               plot(controlll(find(kk),1:trial),controlll(find(jj),1:trial),'o','markeredgecolor','b','markerfacecolor','b','markersize',3)
               hold on
               if k==1, ylabel(strrep(y_factors(j),'_','\_')),set(get(gca,'YLabel'),'Rotation',0.0,'Horizontal','Right','FontWeight','demi'),end
               if j==y_rows, xlabel(strrep(x_factors(k),'_','\_')),set(get(gca,'XLabel'),'FontWeight','demi'),end
               if j==1 & k==ceil(x_cols/2),title(patchlog_name,'FontWeight','bold'),end%(strrep(patchlog_name,'_','\_'),'FontWeight','bold'),end
           end
		end
	end
end

if do_raw_measuring
	%save the full trace from each of the two channels (with azaz), and label it accordingly
	if save_that_sht
        persistent raw_values
		if trial==1, raw_values=[]; end
		
		%padding with zeros if the trial (i.e., the IPI) was shorter than the longest trial in this dataset (hdf stores trials of different lengths, but to put them in a matlab array they must all be the same length)
		if length(ddata(stim_indices(1,1):end,1)) < raw_window(end),   %thus raw_window should be set s.t. the second value covers all the possible range
            ddata=[ddata;zeros((raw_window(end)+1-length(ddata(stim_indices(1,1):end,1))),2)];
		end
		
		for j=1:2
            raw_holding(:,j)=ddata(stim_indices(1,1)+raw_window,j); 
		end
		oldsiz=size(raw_values,2);
		raw_values(1:number_of_IPIs,oldsiz+[1:2])=[azaz' azaz'];
		raw_values(number_of_IPIs+[1:size(raw_holding,1)],oldsiz+[1:2])=raw_holding;
		
		if (trial==number_of_trials) | (save_after_each_trial)
			filename=strcat(common_pathname,'__raw');
			var='raw_values';
			save(filename,var,asky,dub,tbs)
		end
	end
end