function [ sigma, n_end ] = block_half( xx, tt, dt, method )
%
% function [ sigma, n_end ] = block_half( xx, tt, dt, method )
% return parameters of right-half of pulse
% Input:  xx      -- array of counts; xx(1) = value of peak
%         tt      -- array of times;  tt(1) = time of peak
%         dt      -- array of intervals
%         method  -- 1 = fit exponential; 2 = area width
% Output: sigma   -- decay constant
%         n_end   -- index of last block of the pulse
%==============================================================

nx = length( xx );
n_end = 0;
small = 1.e-8;
thresh = 0.05 * xx(1);

for nn = 2:nx

   top_flag = 1;

   if xx( nn ) <= thresh

    % NO DATA!; END OF PULSE
      n_end = nn - 1;
      top_flag = 0;
      break

   elseif xx( nn ) == xx( nn - 1 )

      % FLAT REGION

      if top_flag > 0
         % Keep flat top ...
         n_end = nn;
      else 
         % but not flat bottom
         break
      end

   elseif xx( nn ) < xx( nn - 1 )

      % DECREASING
      top_flag = 0;
      n_end = nn;

   else

      % INCREASING; END OF PULSE
      top_flag = 0;
      break

   end

end
      
if isempty( n_end ) | n_end <= 1

   % No pulse!
   sigma = 0.5*dt(1);
   n_end = 1;

else


   if method == 1

      pulse_amps  = xx(2:n_end)/xx(1);
      pulse_times = tt(2:n_end)-tt(1);
      n_pulse     = length( pulse_amps );

      if n_pulse > 1
         pp = polyfit( pulse_times', log(pulse_amps), 1 );
         log_slope = pp(1);
      else
         log_slope = log( pulse_amps ) / pulse_times;
      end

      if log_slope == 0
         sigma = 1.e99;
      else
         sigma = -1 / log_slope;
      end

   elseif method == 2

      pulse_amps  = xx(1:n_end)/xx(1);
      pulse_times = abs( tt(1:n_end) - tt(1) );
      sigma = sum( pulse_amps * pulse_times ) / sum( pulse_amps );
      if sigma <= 0
         sigma = small;
      end

   end


end