function variance_mean_analysis(events,templates)
%input the templates fitted to the slected events, then vetted for latency (>40), rise time constant (<11), and fit quality (best_resnorm <, say, 5000)
%current plan, wednesday may 10, 5pm -- pull out the events all over again, with more of a tail, then fit with a single exponential for decay, such that long decays can be eliminated

tic
latency_thresh=40 %pick greater than
tau_rise_thresh=11 %pick less than
fit_quality_thresh=5000 %pick less than
tau_decay_thresh_upper=40 %pick less than
tau_decay_thresh_lower=30 %pick greater than

tsize=size(templates),esize=size(events)

picks1=find(templates(:,216)>latency_thresh);   %this index will be different depending on the length of the template and on the number of parameters fitting the template
templates=templates(picks1,:);events=events(picks1,:);
tsize=size(templates),esize=size(events)

picks2=find(templates(:,213)<tau_rise_thresh);
templates=templates(picks2,:);events=events(picks2,:);
tsize=size(templates),esize=size(events)

picks3=find(templates(:,219)<fit_quality_thresh);
templates=templates(picks3,:);events=events(picks3,:);
tsize=size(templates),esize=size(events)

% % figure,plot(templates(:,214),'k.'),title('decays')
picks4=find(templates(:,211)<1);%tau_decay_thresh_upper);
templates=templates(picks4,:);events=events(picks4,:);
tsize=size(templates),esize=size(events)

picks5=find(templates(:,211)>-1);%tau_decay_thresh_lower);
templates=templates(picks5,:);events=events(picks5,:);
tsize=size(templates),esize=size(events)

% figure(20),plot(templates(:,211),'k.'),return
% figure(20),historical_histogram(-1*templates(:,217))

% meanevent=mean(events(:,1:211));
% figure,plot(meanevent)
% meaneventtemplate=iterative_curve_fit_singledecay(meanevent);
% hold on,plot(meaneventtemplate(1:211),'r'),legend(num2str(meaneventtemplate(214)))
% return

%sort the templates according to amplitude (the third-to-last value stored during fitting), and keep the indices of the sort
[sorted_templates_amp sorted_ind_templates_amp]=sort(templates(:,217));

figure,plot(sorted_templates_amp,'k.')

% first_one=input('What is the first event to use?')
% last_one=input('What is the last event to use?')
first_one=60,last_one=140

%select the corresponding raw waveforms, these will be used for v-m analysis
usables=sorted_ind_templates_amp(first_one:last_one);

[minz alignment_points]=min(templates(usables,1:211)');     %finding the peak of the PSC from each template

aligned_template_decays=[];
for vectorize_this=1:length(usables)
	aligned_template_decays=[aligned_template_decays; templates(usables(vectorize_this),alignment_points(vectorize_this)+[0:100])];
end
figure,plot(aligned_template_decays')

%take all the waveforms to be used, take samples from the peak (trough) through n samples later
aligned_event_decays=[];
for vectorize_that=1:length(usables)
	aligned_event_decays=[aligned_event_decays; events(usables(vectorize_that),alignment_points(vectorize_that)+[0:100])];
end
I_macro_vector=mean(aligned_event_decays);
figure, plot(aligned_event_decays'), hold on, plot(I_macro_vector,'r','linewidth',2.5)

% figure,plot(events(usables,1:211)')

baseline_variance_matrix=[];
for need_a_vector=1:length(usables)
    rowww=usables(need_a_vector);
    colll=floor(templates(rowww,216));  %using the latency as found by the template, and going a certain number of samples before that to get baseline
    collls=colll+[-30:-15];
    baseline_variance_matrix=[baseline_variance_matrix;events(rowww,collls)];
end
figure,plot(baseline_variance_matrix')
baseline_variance=mean(var(baseline_variance_matrix))
var_base_var=var(var(baseline_variance_matrix));

%a question:  should the raw waveforms be scaled s.t. their peak amplitude = the mean peak amplitude?  this seems bizarre, yet it is similar to what is done
%a question:  how exactly should the baseline variance be subtracted?
%a question:  should i use the pairwise comparison method of noceti et al. to get a variance for each sample

% var_macro_vector=var(aligned_event_decays);
differences_from_scaled_mean=[];
first_terms_of_mean=mean(I_macro_vector(1:5)');
for okay_then_loopit=1:length(usables)
    decaytrace=aligned_event_decays(okay_then_loopit,:);%figure(10),plot(decaytrace)
    first_terms_of_event=mean(decaytrace(1:5));
    mean_event_decay_scaled=I_macro_vector*(first_terms_of_event/first_terms_of_mean);
    la_difference=decaytrace-mean_event_decay_scaled;
%     figure,plot(decaytrace,'k'),hold on,plot(mean_event_decay_scaled,'r'),legend(num2str(okay_then_loopit))
    differences_from_scaled_mean=[differences_from_scaled_mean;la_difference];
end
var_macro_vector=var(differences_from_scaled_mean);     %should decide whether to bin
keeperz=find(I_macro_vector<-4);I_macro_vector=I_macro_vector(keeperz);var_macro_vector=var_macro_vector(keeperz);
% var_macro_vector=var_macro_vector-baseline_variance;  %baseline variance should be a term in the fitting
figure(9), hold on
plot(-1*I_macro_vector, var_macro_vector,'o','markersize',10,'markeredgecolor','k','markerfacecolor','k')
a=size(I_macro_vector)

%plot the variance versus the mean for each sample -- fit the scatter plot with the inverse parabola sigma^2(t) = <I(t)>i - <I(t)>^2/N -- i is the intitial slope, N is 2I/i 

vm_parabola=vm_curve_fitter(I_macro_vector,var_macro_vector);  %-- consider doing this with brute force
b=size(vm_parabola)
plot(-1*I_macro_vector,vm_parabola,'r')
toc