%
% Program to read in a RTOF file and make some manipulation to the data
% select a .888 file and the corresponding data file we be read
% the program then sorts the data into a nice matrix and rebins as specified in the header
% it also defines a cutoff mass...

function [binned_mass, binned_counts, temperature] = read_rtof()
%----
%adapt here:
binfactor=1;
maxmass=120;
laser_start_frequency=30; %Hz

disp('data will be binned by a factor of: '), disp(binfactor)
disp('maximum mass selected is: '), disp(maxmass)

start_temp=input('Enter start temperature please.')
tempstep=input('Enter temperature increase please (deg K per minute is assumed).')
%------

% first we select the required definition file of the run we want to analyse
[FileName,PathName] = uigetfile('*.888','Select the Re-TOF file describing the run');
fileposition=strcat(PathName, FileName);

% now we have matlab open the file and pull the required info (range and cycles, binwidth):
fileID = fopen(fileposition);
filecontent = textscan(fileID, '%s','delimiter', '\n');
range_string=filecontent{1}{9};
sweeps_string=filecontent{1}{12};
cycle_string=filecontent{1}{15};
bitshift_string=filecontent{1}{19};
range=str2num(range_string(7:end));
cycles=str2num(cycle_string(8:end));
bitshift=str2num(bitshift_string(10));
sweeps=str2num(sweeps_string(10:end));
binwidth=2^bitshift;
fclose(fileID);

%now we have it read the actual data and reformat it to the required matrix
newfileposition=strcat(fileposition(1:(end-3)),'asc');
unsorted_content=importdata(newfileposition);
data=reshape(unsorted_content, range, cycles);

%define a time scale (not the wallclock time but the flight time of the ions): 
time=(1:range)*binwidth;

%this time relates to a mass somewhat according to:
mass= 0.06772 + time*5.41166E-5 + (time.^2)*4.81818E-8 + (time.^3)*-4.7335E-15;
%before alex break the ReTOF values:
%mass= 0.00431 + time*6.95237e-5 + (time.^2)*4.74369e-8 + (time.^3)*3.38808e-15;
%older values:
%mass= -1.05 + time*1.76e-4 + (time.^2)*4.38e-8 + (time.^3)*3.74e-14;

%we record up to a mass of 700 or so. But signal is only up to maybe 150. so I throw away that stuff:
ind_newmass=mass < maxmass;
new_data=data(ind_newmass,:);
new_mass=mass(ind_newmass);

%we have quiet a small binwidth so we want to rebin the whole thing...
[bin_data_row1, bin_mass]=binning_1d(new_mass, new_data(:,1), ceil(size(new_data, 1)/binfactor));
bin_data=zeros(size(bin_data_row1,2), size(new_data,2));
bin_data(:,1)=bin_data_row1;
for kk = 2:size(new_data,2);
	[bin_data(:,kk), bin_mass]=binning_1d(new_mass, new_data(:,kk), ceil(size(new_data, 1)/binfactor));
end;

%output:
binned_counts=bin_data;
binned_mass=bin_mass;
rawdata=data;

% now we make a temperature scale which translates cycle numbers to a temperature scale
%this is how many seconds a cycle takes:
cycletime=sweeps/laser_start_frequency;
time_according_to_cyclenumber=(0:size(new_data,2)-1)*cycletime;
kelvin_per_sec=tempstep/60;
temperature=start_temp+time_according_to_cyclenumber*kelvin_per_sec;

%we can now plot a 2D spectrum (counts versus mass) of a required cycle:
%for now, randomly choose cycle 111: 
%plot(mass, data(:,111))

%or plot all cycles together:
%plot(mass, sum(data,2))

%or plot a range of masses at a range of temperature:
%plot(mass, sum(data,2))

end