function test_dist(filename,units)

plot_s = [1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0];

if nargin < 1
    filename = 'l0_test.txt';
end
if nargin < 2
    units = 'gpt';
    display('assuming gpt units, look out!')
end

mc2 = 510998.946;

switch units
    case 'astra'
        rscale = 1e3;
        pscale = 1e3/mc2;
        tscale = 1e3;
        qscale = 1e3;
        nheaders=0;
        
    case 'gpt'
        rscale = 1e3;
        pscale = 1e3;
        tscale = 1e12;
        qscale = 1e12;
        nheaders=1;
end

data = importdata(filename, ' ', nheaders);


if(isstruct(data))
    data = data.data;
end

x=rscale*data(:,1); % convert [m] to [mm]
y=rscale*data(:,2); % convert [m] to [mm]
z=rscale*data(:,3);
    
px=pscale*data(:,4); % convert [px] to gamma*beta X 1000
py=pscale*data(:,5); % convert [px] to gamma*beta X 1000
pz=pscale*data(:,6); % convert [px] to gamma*beta X 1000

t=tscale*data(:,7); % convert to ps
q=qscale*data(:,8); % convert to pC
nmac = data(:,9);

p = sqrt(px.^2 + py.^2 + pz.^2);
E = sqrt((p*mc2/1e3).^2 + (mc2)^2);

Qbunch = sum(q.*nmac);

avgx = mean(x);
stdx = std(x,1);

avgy = mean(y);
stdy = std(y,1);

avgz = mean(z);
stdz = std(z,1);

stdpx = std(px,1);
stdpy = std(py,1);
stdpz = std(pz,1);

avgt = mean(t);
stdt = std(t,1);

avgE = mean(E);
stdE = std(E);
dE = stdE/avgE;
dP = std(p)/mean(p);

enx = stdx*stdpx;  % en = stdx * (c*stdpx) / mc^2 if [stdpx] =eV/c
eny = stdy*stdpy;


GB=

kTx = 0.511e6*(stdpx/1e3)^2;
kTy = 0.511e6*(stdpy/1e3)^2;
%kT = (kTx+kTy)/2;

%pinhole_size_x = abs(min(x))+max(x);
%pinhole_size_y = abs(min(y))+max(y);
%pinhole_size = (pinhole_size_x+pinhole_size_y)/2;

if plot_s(1)==1
    figure
    plot(x,y,'*')
    xlabel('x (mm)')
    ylabel('y (mm)')
    title({['<x> = ' num2str(avgx) ' (mm), ' '\sigma_x = ' num2str(stdx) ' (mm)'],['<y> = ' num2str(avgy) ' (mm), \sigma_y = ' num2str(stdy) ' (mm)']})
    axis equal
    
    f=1/2;
    nbins = floor(f*sqrt(length(x)));
    
    %r = sqrt();
    figure
    rho = hist3([x y],[nbins,nbins]);
    rho=rho/max(max(rho));
    imagesc(rho);
    colorbar
    %figure
    %radial_density_histogram_plot(x, y, ones(size(x)), 500)
end

if plot_s(2)==1
    
    figure
    ft=1;
    nbins=floor(ft*sqrt(length(t)));
    
    [rho_t, tcs]=hist(t,nbins);
    
    dt = tcs(2)-tcs(1);
    Qtot = sum(rho_t*mean(q));
    
    plot(tcs,rho_t*mean(q)/dt)
    xlabel('emission time (ps)')
    ylabel('I (A)')
    title(['Qbunch = ' num2str(abs(sum(q))) ' (pC), \sigma_t = ' num2str(stdt) ' (ps)'])
end

if plot_s(3)==1

  figure
  subplot(1,2,1)
  plot(x,px*mc2/1e3,'*')
  xlabel('x (mm)')
  ylabel('p_x (eV/c)')
  title(['\epsilon_{n,x} = ' num2str(enx) ' (mm-mrad), kT = ' num2str(kTx*1e3) ' (meV)'])
  
  subplot(1,2,2)
  plot(y,py*mc2/1e3,'*')
  std(py*mc2/1e3)
  xlabel('y (mm)')
  ylabel('p_y (eV/c)')
  title(['\epsilon_{n,y} = ' num2str(eny) ' (mm-mrad), kT = ' num2str(kTy*1e3) ' (meV)'])

end

if plot_s(4)
    figure
    plot(z,pz*mc2/1e3,'*')
    xlabel('z (mm)')
    ylabel('p_z (eV/c)')
end

if plot_s(5)
    %figure
    %plot(t,sqrt(x.^2+y.^2),'*')
    %xlabel('t (ps)')
    %ylabel('r (mm)')
    
    figure
    plot(t,z,'*')
    xlabel('t (ps)');
    ylabel('z (mm)');
    
    figure
    plot(t,x,'*')
    xlabel('t (ps)')
    ylabel('x (mm)')
    
    figure
    plot(t,y,'*')
    xlabel('t (ps)')
    ylabel('y (mm)')
    
    figure
    plot(t,E,'*')
end

display(['Qbunch = ' num2str(Qbunch) ' (pC)'])
%display(['pinhole size = ' num2str(pinhole_size) ' (mm)'])
display(['avgx = ' num2str(avgx) ,' (mm), avgy = ' num2str(avgy) ' (mm), avgz = ' num2str(avgz) ' (mm)'])
display(['stdx = ' num2str(stdx) ,' (mm), stdy = ' num2str(stdy) ' (mm), stdz = ' num2str(stdz) ' (mm)'])
display(['avgt = ' num2str(avgt) ' (ps), stdt = ' num2str(stdt) ' (ps)'])
display(['avgE = ' num2str(avgE/1e6) ' (MeV), dE = ' num2str(dE) ' ()'])
display(['norm. emittance x = ' num2str(enx) ' (mm-mrad), kT = ' num2str(1e3*kTx) ' (meV)'])
display(['norm. emittance y = ' num2str(eny) ' (mm-mrad), kT = ' num2str(1e3*kTy) ' (meV)'])

end




