function tumorBeamletDVH(lin245cm,lin100cm,lin0625cm,lin025cm)

addpath('C:\Documents and Settings\Steve\My Documents\Research\matlab')
addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqUtil')
addpath('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input')

Ndvhbins = 1000;                                                            %number of dvh bins

%open prescription in Matlab
%prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\HN003linPrescVol\PTV70DosePlus.bin');
%prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\powerLaw\linPrescJake039cm\tumorDosePlus.bin');
prescTumorLin = open_presc('C:\Documents and Settings\Steve\My Documents\Research\optNew\linlsqOpt\input\HN003powerLaw\normPower1\PTV70DosePlus.bin');


%dim = 1;                                                                    %2D
dim = 30;                                                                   %3D

fid = fopen(lin245cm,'rb');                                           %open plans in Matlab
planLin0625cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(lin100cm,'rb');                                           %open plans in Matlab
planRoot0625cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(lin0625cm,'rb');                                           %open plans in Matlab
planQuad0625cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

fid = fopen(lin025cm,'rb');                                           %open plans in Matlab
planGomp0625cm = reshape(fread(fid,'float'),72,72,dim);
fclose(fid);

tumorMaskLin = single(~~(full3D(prescTumorLin)));                                 %define tissue masks
tumorMaskRoot = single(~~(full3D(prescTumorLin)));
tumorMaskQuad = single(~~(full3D(prescTumorLin)));
tumorMaskGomp = single(~~(full3D(prescTumorLin)));


tumorDoseLin0625cm = planLin0625cm.*tumorMaskLin;                                                %dose distributions in each ROI
tumorDoseRoot0625cm = planRoot0625cm.*tumorMaskRoot;
tumorDoseQuad0625cm = planQuad0625cm.*tumorMaskQuad;
tumorDoseGomp0625cm = planGomp0625cm.*tumorMaskGomp;


%convert 3D matrices to sparse matrices
tumorDoseSparseLin0625cm = sparse3D(tumorDoseLin0625cm);
tumorDoseSparseRoot0625cm = sparse3D(tumorDoseRoot0625cm);
tumorDoseSparseQuad0625cm = sparse3D(tumorDoseQuad0625cm);
tumorDoseSparseGomp0625cm = sparse3D(tumorDoseGomp0625cm);

tumorNormDoseLin0625cm = ((tumorDoseSparseLin0625cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseLin0625cm = mean(prescTumorLin.non_zero_values);
tumorEdoseLin0625cm = tumorNormDoseLin0625cm.*tumorAverageDoseLin0625cm;
tumorEdoseLin0625cm = tumorNormDoseLin0625cm;

tumorNormDoseRoot0625cm = ((tumorDoseSparseRoot0625cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseRoot0625cm = mean(prescTumorLin.non_zero_values);
tumorEdoseRoot0625cm = tumorNormDoseRoot0625cm.*tumorAverageDoseRoot0625cm;
tumorEdoseRoot0625cm = tumorNormDoseRoot0625cm;

tumorNormDoseQuad0625cm = ((tumorDoseSparseQuad0625cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseQuad0625cm = mean(prescTumorLin.non_zero_values);
tumorEdoseQuad0625cm = tumorNormDoseQuad0625cm.*tumorAverageDoseQuad0625cm;
tumorEdoseQuad0625cm = tumorNormDoseQuad0625cm;

tumorNormDoseGomp0625cm = ((tumorDoseSparseGomp0625cm.non_zero_values)./(prescTumorLin.non_zero_values));  %calculate equivalent uniform dose
tumorAverageDoseGomp0625cm = mean(prescTumorLin.non_zero_values);
tumorEdoseGomp0625cm = tumorNormDoseGomp0625cm.*tumorAverageDoseGomp0625cm;
tumorEdoseGomp0625cm = tumorNormDoseGomp0625cm;

dMax = 1.1*max(tumorEdoseRoot0625cm);
dvhbins = [0:Ndvhbins-1]*dMax/Ndvhbins;

dHistTumorLin0625cm = histc(tumorEdoseLin0625cm,dvhbins);
dHistTumorRoot0625cm = histc(tumorEdoseRoot0625cm,dvhbins);
dHistTumorQuad0625cm = histc(tumorEdoseQuad0625cm,dvhbins);
dHistTumorGomp0625cm = histc(tumorEdoseGomp0625cm,dvhbins);

dHistTumorNormLin0625cm = 100*dHistTumorLin0625cm./(double(tumorDoseSparseLin0625cm.Nind));
dHistTumorNormRoot0625cm = 100*dHistTumorRoot0625cm./(double(tumorDoseSparseRoot0625cm.Nind));
dHistTumorNormQuad0625cm = 100*dHistTumorQuad0625cm./(double(tumorDoseSparseQuad0625cm.Nind));
dHistTumorNormGomp0625cm = 100*dHistTumorGomp0625cm./(double(tumorDoseSparseGomp0625cm.Nind));

DVHtumorLin0625cm = flipud(cumsum(flipud(dHistTumorNormLin0625cm)));
DVHtumorRoot0625cm = flipud(cumsum(flipud(dHistTumorNormRoot0625cm)));
DVHtumorQuad0625cm = flipud(cumsum(flipud(dHistTumorNormQuad0625cm)));
DVHtumorGomp0625cm = flipud(cumsum(flipud(dHistTumorNormGomp0625cm)));

% %rescale 1.00 cm due to misregistration
DVHtumorRoot0625cm = DVHtumorRoot0625cm./(dvhbins.')-1.8719;

tumorfigure(dvhbins,DVHtumorLin0625cm,DVHtumorRoot0625cm,DVHtumorQuad0625cm,DVHtumorGomp0625cm)

end