% function accidental(f1, f2, gate1, gate2)
%    f1: rate of channel-1
%    f2: rate of channel-2
%    gate1: signal-width that goes into coincidence unit (channel-1)
%    gate2: signal-width that goes into coincidence unit (channel-2)
function accidental(f1, f2, gate1, gate2)


% last trigger time for each channel (left-right)
t1=0;
t2=0;


% next event time
% Next e-time probablity is determined by Poisson distribution -> P =  f*exp(-t*f)
t1_next = t1 - log(rand(1))/f1;
t2_next = t2 - log(rand(1))/f2;


% current state of each side (logic gates)
s1=0;
s2=0;


% delta-t array (this will contain (t2-t1) triggers)
dT=zeros(1,1000000);
n=0; % size of array
n_pre = 0;


% current time
t=0;


% historam of (t2-t1) is limited by the following
t_min = -gate2*1.1;
t_max = gate1*1.1;


% make an infinite loop
% in  this loop determine what happens in the "future"
while 1

  % Draw the histogram every given coincidence event
  if (mod(n,10)==0 & n_pre~=n)
    n_pre = n;
    hist(dT(1:n)*1e9, (t_min:(t_max-t_min)/50:t_max)*1e9);
    drawnow
  end

  
  % there are four possible events in the future
  %   1) The end of t1 pulse (gate1 width) : t1+gate1
  %   2) The end of t2 pulse (gate2 width) : t2+gate2
  %   3) Next t1 pulse                  : t1_next
  %   4) Next t2 pulse                  : t2_next

  % We'll choose the earliest one, and change the state accordingly
  % note that 1) and 2) can happen only if the state of that channel is 
  % logic-1.

  % Once we determine the earliest case, we take the appropiate
  % actions and then 'continue' to the next step (loop)
  

  % 1) is it end of t1 ?
  if (s1==1) & (s2==0 | (t1+gate1) < (t2+gate2)) & ((t1+gate1) < t1_next) & ((t1+gate1)<t2_next)
    % change state-s1 to logic-0, and set the current time
    s1 = 0;
    t = t1+gate1;
    continue;
  end

  
  % 2) is it end of t2 ?
  if (s2==1) & (s1==0 | (t2+gate2) < (t1+gate1)) & ((t2+gate2) < t1_next) & ((t2+gate2)<t2_next)
    % change state-s2 to logic-0, and set the current time
    s2 = 0;
    t = t2+gate2;
    continue;
  end

  
  % 3) is it time for next t1
  if (s1==0 | t1_next < (t1+gate1)) & (s2==0 | t1_next < (t2+gate2)) & (t1_next < t2_next)
    % Reset the t1 time, and determine the next trigger of this
    % channel according to the Poisson distribution.
    % Make if state changes from 0->1 and the s2 is already logic-1
    % it is a coincidence trigger -> calculate (t2-t1) and store it
    % into the array. (It deselects the pileup. Move s1=1 before
    % "if", and remove (s1==0), to ignore the pileup.)
    t1 = t1_next;
    t1_next = t1 - log(rand(1))/f1;
    t = t1;
    if (s1==0) & (s2==1)
      n = n+1;
      dT(n) = t2-t1;
    end
    s1 = 1;
    continue;
  end

  
  % 4) is it time for next t2
  if (s1==0 | t2_next < (t1+gate1)) & (s2==0 | t2_next < (t2+gate2)) & (t2_next < t1_next)
    % Reset the t2 time, and determine the next trigger of this
    % channel according to the Poisson distribution.
    % Make if state changes from 0->1 and the s1 is already logic-1
    % it is a coincidence trigger -> calculate (t2-t1) and store it
    % into the array. (It deselects the pileup. Move s2=1 before
    % "if", and remove (s2==0), to ignore the pileup.)
    t2 = t2_next;
    t2_next = t2 - log(rand(1))/f2;
    t = t2;
    if (s1==1) & (s2==0)
      n = n+1;
      dT(n) = t2-t1;
    end
    s2 = 1;
    continue;
  end

  
  % It comes here if there is a bug in above code
  disp('Error: It should not come to this point');
  [t1_next  t2_next  s1*(t1+gate1)  s2*(t2+gate2)]
  t1
  t2
  t1_next
  t2_next
  s1
  s2
  t
  whos
  return

end

