clear all

[filename_,PathName]=uigetfile('*.mp','Multiselect','on');
%[FileName,PathName] = uigetfile('*.888','Select the Re-TOF file describing the run'); %Select the TOF file you want to load
if isnumeric(filename_)
    disp('File selection aborted by user')
    return
end

if iscell(filename_)
   FileName=filename_; 
else
   FileName=cellstr(filename_);%in case of one element filename_ has to be converted to a cell
end
for i=1:length(FileName)

close all
max_mass=200;
laser_start_frequency=30;
start_temp=0;
tempstep=1;
fileposition=strcat(PathName, FileName{1,i});
a=importdata(fileposition,' ',54);
range_string=a.textdata{2};
sweeps_string=a.textdata{7};
cycle_string=a.textdata{12};
bitshift_string=a.textdata{43};
range=str2num(range_string(7:end));
cycles=str2num(cycle_string(8:end));
bitshift=str2num(bitshift_string(10));
sweeps=str2num(sweeps_string(10:end));
% fclose(fileID);
% newfileposition=strcat(fileposition(1:(end-3)),'asc');
%binfactor=2;
binwidth=0.8*2^bitshift;

time=(0:range-1)*binwidth;
a=importdata(fileposition,' ',54);
data = reshape(a.data,range,cycles);
%data(:,105)=data(:,104);
designation=FileName{1,i}(1:end-4);
%image(data);
 %d = b*10;
 %image(d);
 %3D data
%save('CO_H2Orich_HigDos_10K_10_49eV.txt','data','-ascii');
c=sum(data,2);
c=smooth(c,5);
%figure, plot(c);
%mass= 1.007*(0.05746885+ time*5.668314e-5+(time.^2)*4.758394e-8);
% mass_cal_factor=1.006; %% for blank 10.49
% mass_cal_factor=1.004; %% for normal irradiation 10.49
% mass_cal_factor=1.003; %% for CD3OH:CD3CDO irradiation 10.49
% mass_cal_factor=1.002; %% for CD3OH:CD3CDO irradiation 9.70
% mass_cal_factor=1.004; %% for CH3OH:CD3CHO blank 9.70
% mass_cal_factor=1.003; %% for CH3OH:CD3CHO irradiation 9.70
% mass_cal_factor=1.004; %% for CH3OH:CH3CHO irradiation 9.70
% mass_cal_factor=1.002; %% for CH3OH:CH3CHO blank 9.70
% mass_cal_factor=1.004; %% for CH3OH:CH3CHO irrad 9.20
% mass_cal_factor=1.004; %% for CD3OH:CD3CDO irrad 9.20
% mass_cal_factor=1.004; %% for CH3OH:CD3CHO irrad 9.20
% mass_cal_factor=1.0055; %% for CH3OH:CD3CHO irrad 8.25
% mass_cal_factor=1.006; %% for CH3OH:CD3CHO irrad 8.64
% mass_cal_factor=1.004; %% for CH3OH:CD3CHO irrad 
mass_cal_factor=1.0025; %% for H2S:CO irrad 10.82
% mass_cal_factor=1.006; %% for H2S:13CO irrad 10.82

 mass= mass_cal_factor*(0.05746885+ time*5.668314e-5+(time.^2)*4.758394e-8);

% % Mass assignment parameters
% m1 = 80; 
% bin1 = 12605; 
% m2 = 52;
% bin2 = 10122;

% Mass assignment parameters
% Mass assignment parameters


% % Mass assignment parameters
% m1 = 2; 
% bin1 = 1823; 
% m2 = 96;
% bin2 = 13843;
% % Mass assignment parameters
% m1 = 1; 
% bin1 = 1232; 
% m2 = 60;
% bin2 = 10895;
% Mass assignment parameters
m1 = 32; 
bin1 = 7889; 
m2 = 98;
bin2 = 13957;


% Mass assignment
t1 = bin1 * binwidth;
t2 = bin2 * binwidth;
A = (m1+m2)/(t1-t2)^2-2*sqrt(m1*m2)/abs(t1-t2)^2;
B = sqrt(m2/A)-t2;
mass = A*(time+B).^2;

plot(mass,c)
findpeaks(c,'MinPeakDistance',35,'MinPeakProminence',5)
[~,Peaks]=findpeaks(c,'MinPeakDistance',35,'MinPeakProminence',5);
%figure
%plot(mass,c)
%masses(1,2)=input('Enter lower mass');
%masses(1,1)=input('Enter channel for lower mass');
%masses(2,2)=input('Enter intermediate mass');
%masses(2,1)=input('Enter channel for intermediate mass');
%masses(3,2)=input('Enter higher mass');
%masses(3,1)=input('Enter channel for higher mass');
[~,Reference]=max(c(1500:end));
mass_ref=mass(Reference+1499);
mass_ref=round(mass_ref);
%disp 'I am guessing your maximum peak is at mass', disp(mass_ref), 
%mass_corr=input('is that correct? y/n','s');
%if mass_corr=='y'
%mass_correction=mass_ref;
%else if mass_corr=='n'
%        mass_correction=input('please enter correct mass');
%    else
%        disp 'error, invalid input'   
%        return
%    end
%end
Peaks_mass=mass(Peaks)';
Peaks_mass=round(Peaks_mass);
%Peaks_mass(1)=2;
TIME=time(Peaks);
%fo = fitoptions('StartPoint',[2.555e-27,5.042e-22,3.824e-17,1.372e-12,2.552e-08,0.0001438,0.2733]);
%f=fit(TIME',Peaks_mass,'poly8')

figure
plot(TIME,Peaks_mass,'o')
hold on
plot(mass,c)
min_mass=1;
hold off
ind_newmass=mass < max_mass;
new_data=data(ind_newmass,:);
new_mass=mass(ind_newmass);
counts=new_data;
mass=new_mass;

figure
plot(mass,sum(counts,2))
% 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;

% 
% %On/Off switches
% Waterfall_Bin = 0; % Binning along mass axis for waterfall plots, set to 0 to turn off
% VUV_Correction = 1; % set to 1 to generate VUV normalization constant for each point in TPD profile
% Generate_EPS = 1; % set to 1 to generate EPS file containing TPD profiles
% 
% smoothing = 1; % # of points used for Boxcar averaging of TPD profiles
% VUV_TStart = 57; % Temperature Tek Acquire was turned on
% VUV_TEnd = 320; % Temperature Tek Acquire was turned off, usually 320 K
% 
% % VUV power correction calculations
% if VUV_Correction == 1
% [FileName2,PathName2] = uigetfile('*.txt');
% VUV = importdata(strcat(PathName2,FileName2));
% VUV = VUV(:,2);
% T_VUV = linspace(VUV_TStart,VUV_TEnd,length(VUV));
% T_VUV = interp1(temp,temp,T_VUV,'nearest');
% Export_VUV = zeros(length(temp),1);
% for i = (1:length(temp))
%   TIndex = T_VUV == temp(i);
%   VUV_Intensity = mean(VUV(TIndex,:));
%   Export_VUV(i) = VUV_Intensity;
% end
% Export_VUV = Export_VUV/1000;
% Export_VUV(isnan(Export_VUV)) = 1;
% else
% Export_VUV = ones(length(temp),1);
% end


number_specs=max_mass-min_mass;
mass_ind = (mass > min_mass & mass < max_mass);
toplot=smooth(sum(counts,2),10);
%semilogy (mass(mass_ind), toplot(mass_ind))  ;
peakwidth=0.8; % changes the range around which mass center the passes are integrated. (should be decreased if non-integer m/c is expected)
min_set=10; % sets a minimum threshold to plot 
cycletime=sweeps/laser_start_frequency;
time_according_to_cyclenumber=(0:size(new_data,2)-1)*cycletime;
kelvin_per_sec=tempstep/60;
temp=start_temp+time_according_to_cyclenumber*kelvin_per_sec;

figure
smooth_expdata=[];

for kk = min_mass:max_mass;
	mass_ind = (mass > kk - peakwidth/2 & mass < kk + peakwidth/2);
    smoth_data=smooth( sum(counts(mass_ind,:),1) ,1);
    smooth_expdata=[smooth_expdata,smoth_data];
 	plot (temp, smooth( sum(counts(mass_ind,:),1) ,1))  ;
	%hold on
	%plot (temp, smooth( sum(counts2(mass_ind,:),1) ,5), 'r')  ;
	
    xlabel('T (K)');
		ylabel('counts');
		title(strcat({'mass: '},{num2str(kk)},{' amu'}, {' '}, {designation}),'Interpreter','none');
		y_axis_set=max(smooth( sum(counts(mass_ind,1:end),1) ,5));
		if y_axis_set > min_set ;
			y_max=y_axis_set+0.1*y_axis_set; %set ymax of plot
		else
			y_max=min_set; % if nothing is there, set a ymax anyway to not be confused with noise
        end
		axis([1 320 0 y_max]);
		% set(gca,'ytick',[])
		% set(gca,'yticklabel',[])
		set(gcf, 'PaperUnits', 'inches');
		set(gcf, 'PaperType', 'usletter');
		set(gcf, 'PaperPosition', [0.5 0.5 8.5 4.5]);
		set(gca,'FontSize',20,'LineWidth',1.25)
	if kk == min_mass
		print ('-dpsc', strcat(designation, '_single_masses_over_temp.eps'));
	else
		print ('-dpsc', strcat(designation, '_single_masses_over_temp.eps'), '-append');
	end
	%hold off
end
end
%eps2pdf(strcat(designation, '_single_masses_over_temp.eps'))