function out = lagHist(data1,data2,range,max)
% スパイク時刻差ヒストグラム
% [負 0 正] = lagHist(データ1,データ2,試行範囲,最大時間差)
%
% データ1のスパイク {t1_i}とデータ2のスパイク {t2_i}との時間差 
% l_ij = t1_i - t2_i のヒストグラムをデータ精度と同じビンで返す関数
%
% out(:,1) : l= -1,-2, ..., -Max のヒストグラム
% out(1,2) : l=0 のヒストグラム
% out(:,3) : l= 1,2, ..., Max のヒストグラム
%


trial=size(range,1);		%%% 試行数
len=range(1,2)-range(1,1);	%%% 切り出す長さ

lhist=zeros(2*max+1,1);
num=zeros(1,2);

for i=1:trial
  train1=data1(range(i,1)<data1 & data1<=range(i,1)+len);
  train2=data2(range(i,1)<data2 & data2<=range(i,1)+len);
			%%% 試行ごとに切り出す．
  num = num + [length(train1) length(train2)];

  for j1=1:length(train1)
     for j2=1:length(train2)
       lag=train1(j1)-train2(j2);
	   %%% data2のスパイク時刻を基準とした、data1のスパイクの相対時刻列

	if(-max<=lag & lag<=max)
           lhist(max+lag+1) = lhist(max+lag+1)+1/(len-abs(lag));
        end
     end
  end
end

lhist = lhist*trial*len*len/num(1)/num(2);

out= zeros(max,3);
out(:,1)= lhist(max:-1:1) ;
out(1,2)= lhist(max+1);
out(:,3)= lhist( (max+2):(2*max+1) );