function best_intervals=montecarlo_noisy_events(runtimes)
%simulates (that is, is a simulacrum of) recording noise and variable synaptic release
%we're interested in how many are released
%so we ask:  do the notches in the rising and falling phase suggest multiple release events, or just chance alignments of dropping noise + rising signal

%how many trials
%what is the outcome to be measured, e.g., how often do you see n consecutive events where the mean is linear and horizontal? (i.e., i guess, where the points will be well...
%fit (how well?) with a line with a slope between .x and -.x)

% variance reduction, antithetic variates, control variates, importance sampling, stratified sampling

% compare single peak versus multipeak to see which is more likely

close all, tic
amprange1=[4 8]
tauriserange1=[15 20]
taudecrange1=[50 100]

second_release=0
amprange2=[4 8]
tauriserange2=[15 20]
taudecrange2=[50 100]
tem2delayrange=[1 50]

noiserange=[2 6]

begin_taper=.8
taperlength=50

% runtimes=100


time_vector=[1:200]; 
best_intervals=[]
for runna=1:runtimes
    %making the template
	amplitude_scaler1=randint(1,1,amprange1);
	tau_rise1=randint(1,1,tauriserange1);
	tau_decay1=randint(1,1,taudecrange1);
	template=-1*[amplitude_scaler1*(1-exp(-time_vector/tau_rise1)).*exp(-time_vector/tau_decay1)];
	plot(template,'k'),hold on
	%tapering the beginning and the end
	a=repmat(begin_taper,1,taperlength); b=[(length(a)+1):-1:2]; c=a.^b; d=c*template(1); template=[d template];
	ratte=(template(end)/template(end-1));
	a=repmat(ratte,1,taperlength); b=[2:(length(a)+1)]; c=a.^b; d=c*template(end); template=[template d];
	plot(template,'r'),hold on
	
	%if indicated, making the second template
	if second_release
		%tem2=.5*template;
		amplitude_scaler2=randint(1,1,amprange2);
		tau_rise2=randint(1,1,tauriserange2);
		tau_decay2=randint(1,1,taudecrange2);
		template2=-1*[amplitude_scaler2*(1-exp(-time_vector/tau_rise2)).*exp(-time_vector/tau_decay2)];
		%tapering the beginning and the end
		a=repmat(begin_taper,1,taperlength); b=[(length(a)+1):-1:2]; c=a.^b; d=c*template2(1); template2=[d template2];
		ratte=(template2(end)/template2(end-1));
		a=repmat(ratte,1,taperlength); b=[2:(length(a)+1)]; c=a.^b; d=c*template2(end); template2=[template2 d];
		tem2delay=zeros(1,randint(1,1,tem2delayrange));
		template2=[tem2delay template2];
        %adding the two templates
        template=template+template2(1:length(template));        %because here we are interested in the rise phase, we don't worry about decay for now
%         plot(template,'g'),hold on
	end
	
	[y,min_ind]=min(template);
	
	%make the noise
	noise=rand(1,length(template));noise=noise*randint(1,1,noiserange);noise=noise-mean(noise);
	%add the noise to the templates (note that this is *after* the templates are added together, which is likely more biological)
	ntemplate=template+noise;
	plot(ntemplate,'b')
	
	%readout:  find the two points during the rising phase st the later point is greater than the former point, separated by the greatest interval 
	%do the same thing on the experimental traces
	
	mm=1;   %two different judging methods -- choose one
% 	for j=taperlength:min_ind
% %         gg=ntemplate(j)
%         afterz=ntemplate((j+1):min_ind);
% %         greaterz=afterz>ntemplate(j)
% %         bin_greaterz=find(afterz>ntemplate(j))
%         long_bin_greaterz=max(find(afterz>ntemplate(j)));
%         if long_bin_greaterz>mm,mm=long_bin_greaterz;end
% 	end
    
    for j=taperlength:min_ind-1
        for k=min_ind:-1:j
            if mean(ntemplate(j:k))>=ntemplate(j)
                longg=k-j;
                break
            end
        end
        if longg>mm,mm=longg;end
    end
    
	best_intervals=[best_intervals mm];
end
toc