function par_mat = block_pulses( t_seg_1, t_seg_2, xx_vec, method )
% function par_mat = block_pulses( t_seg_1, t_seg_2, xx_vec, method )
% Do crude deconvolution of output from make_segments
% Input: t_seg_1 --  array consisting of the times at left end of segments
% Input: t_seg_2 --  array consisting of the times at right end of segments
%        xx_vec  --  array consisting of the mean rates in the segments
%        method  --  method for "block_half"
% Output: par_mat( num_pulses, : ) -- matrix of pulse parameters
%                                  -- first index = pulse number
%                                  -- second index: amp, t_max, sig_left, sig_right, 
%                                                   nn_peak, nn_left, nn_right
%
%------------------------------------------------------------------------------------
new_method = 0;

if new_method > 0

else

   max_xx = max( xx_vec );
   n_x = length( xx_vec );

   times_use = ( t_seg_1 + t_seg_2 ) / 2;
   dt_use    = ( t_seg_2 - t_seg_1 );

   num_pulses = 0;
   
   [ xx_amp, nn_max ] = max( xx_vec );
   par_mat = [];
   thresh = 0.05 * xx_amp; % don't include weak pulses 
   % these are probably background fluct's anyway

   while xx_amp > thresh

      num_pulses = num_pulses + 1;
      tt_max = times_use( nn_max );

      %---------------------------------------
      % Fit right-hand part of pulse
      %---------------------------------------

      if nn_max >= n_x 

         % Pulse peaks at the last point
         sigma_right = dt_use( n_x ); % presumes sharp decline
         nn_right    = nn_max;

      else

         xx =    xx_vec( nn_max: n_x );
         tt = times_use( nn_max: n_x );
         dt = dt_use( nn_max: n_x );

         [ sigma_right, nn_end ] = block_half( xx, tt, dt, method );
         nn_right = nn_max + nn_end - 1; % 1 --> 2?

      end
      if nn_right < nn_max;nn_right = nn_max;end

      %---------------------------------------
      % Fit left-hand part of pulse
      %---------------------------------------

      if nn_max <= 1 

         % Pulse peaks at the first point
         sigma_left = dt_use(1); % presumes sharp decline
         nn_left    = nn_max;

      else

         xx =    xx_vec( nn_max:-1:1 );
         tt = times_use( nn_max:-1:1 );
         dt = dt_use( nn_max:-1:1 );

         [ sigma_left, nn_end ] = block_half( xx, tt, dt, method );
         sigma_left = abs( sigma_left );

         nn_left = nn_max - nn_end + 1; % 2 instead of 1, so that ...

      end

      %-------------------------
      % Compute area of pulse  %
      %-------------------------

      area_vec = xx_vec( nn_left:nn_right ) * ...
              ( t_seg_2( nn_left:nn_right ) - t_seg_1( nn_left:nn_right ) );
      pulse_area = sum( area_vec );

      if nn_left > nn_max;nn_left = nn_max;end

      par_mat( num_pulses, 1 ) = xx_amp;
      par_mat( num_pulses, 2 ) = tt_max;
      par_mat( num_pulses, 3 ) = sigma_left;
      par_mat( num_pulses, 4 ) = sigma_right;
      par_mat( num_pulses, 5 ) = nn_max;
      par_mat( num_pulses, 6 ) = nn_left;
      par_mat( num_pulses, 7 ) = nn_right;
      par_mat( num_pulses, 8 ) = pulse_area;

      xx_vec ( nn_max ) = -1;
      xx_vec( nn_left:nn_right ) = - ones( size( nn_left:nn_right ) );
      [ xx_amp, nn_max ] = max( xx_vec );

   end

end