function [d, thresholds, true1s, false1s, true2s, false2s, pratio] = discriminability(spikes1, spikes2, tmax, T, cmax)
[c1, b1] = counthist(spikes1, 0.0, tmax, T, cmax);
[c2, b2] = counthist(spikes2, 0.0, tmax, T, cmax);
thresholds = 0:cmax;
true1s = zeros(length(thresholds), 1);
true2s = zeros(length(thresholds), 1);
false1s = zeros(length(thresholds), 1);
false2s = zeros(length(thresholds), 1);
for k = 1:length(thresholds)
    th = thresholds(k);
    t1 = sum(c1(b1<=th));
    f1 = sum(c1(b1>th));
    t2 = sum(c2(b2>=th));
    f2 = sum(c2(b2<th));
    true1s(k) = t1;
    true2s(k) = t2;
    false1s(k) = f1;
    false2s(k) = f2;
end
%pratio = (true1s + true2s)./(false1s+false2s);
pratio = (true1s + true2s)/2;
d = max(pratio);
end