%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: 
% -------
%	Dr. Behrouz Eslami-Mossallam
%	Delft University of Technology
%	2629HZ Delft
%
% This script constructs a dwell time distribution of tanslocating polymerases, optimized for data acquired with Magnetic or Optical Tweezers. 
%
% Input variables:
% ----------------
% x_array = the array of enzyme position, converted to nucleotides [nt] or base pairs [bp].
% t_array = the array of timepoints; has to have the same length as x_array.
% Ts =  width of the filtering window in [nt] or [bp]; it is the number of points (odd number) over which the trace is being averaged.
% Nw =  width of the dwell time window. Since that x_array is given in the unit of [nt] or [bp], best choose an integer. 
% Ti = for histogram; Ti is the location of the first bin edge. 
% bins_per_decade = for histogram; defines number of bins per each decade. 
% number_of_decades = for histogram; defines number of decades spanned by the histogram. 
% correction = binary variable (0 or 1). 1: bin width is corrected to be a multiple of the temporal resolution (i.e. time between two consecutive data points).
% nboot = number of bootstrapped samples (integer). Bootstrapping is used to estimate the errorbars in the dwell time distribution.
%
% Example:
% --------
% Example.m = example script that computes and plots the dwell time distribution for the accompanied trace data example (DataExample.mat).
%
% Citation for use of the code:
% -----------------------------
% Richard Janissen, Behrouz Eslami-Mossallam, Irina Artsimovitch*, Martin Depken*, Nynke H. Dekker*
% High-throughput single-molecule experiments reveal heterogeneity, state-switching, and three interconnected pause states in bacterial transcription.
% Cell Reports    (2022)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [bar_pos, DWT_hist, Err] = DWT_analysis(x_array, t_array, Ts, Nw, Ti,bins_per_decade,number_of_decades,correction, nboot)

	x_array = double(x_array);

	t0 = mean(diff(t_array)); % time resolution

	xs_array = sgolayfilt(x_array,0,Ts);

	[Tc_array]=find_crossing_time(t_array, xs_array,Nw);

	DWT_array = diff(Tc_array);

	[bar_pos, bins]=make_bins(Ti,bins_per_decade,number_of_decades,t0,correction);

	[DWT_hist] = DwellTimeHist_v3(DWT_array, t0, bins);

	bootstat = bootstrp(nboot,@ (DWT_array) DwellTimeHist_v3(DWT_array, t0, bins),DWT_array);

	Err=zeros([1, length(DWT_hist)]);

	for i=1:length(DWT_hist)

		Err(:,i)=std(bootstat(:,i));

	end


