WiscPlan_v2/WiscPlanPhotonkV125/matlab_frontend/helicalDosecalcSetup7_fullRO.m

%%

% This version was modified by Peter Ferjancic to accept .nrrd file format
% for target. It was also rollbacked on some parameters back to the
% clinical system settings.

% This file has been modified by Surendra Prajapati especially to run
% WiscPlan for KV beams. Works for running locally as well as in Server.
%
% RTF 11/4/05
% Sets up the run files needed for running a Condor convolution
% superposition dose calculation with photon beams.

% Versions of the convolution code that run in both Windows XP, Condor, and
% directly through Matlab can all be created here.  The only real
% differences between the files are in the access conventions, which just%
% This has to do with switching forward slashes with backslashes.

% Surendra edits: num_batches, SAD, pitch, xpmin/max, ypmin/max, Mxp, Nphi 
% Also edit: 
% Kernel should always be named Kernels.mat 

function num_batches = helicalDosecalcSetup7_fullRO(patient_dir, OptGoals, beamlet_dir)
% -- INPUT:
% patient_dir: specify where kernel and [geometry] files are located

iso = [0 0 0]; % Point about which gantry rotations begin
SAD = 85;  % source-axis distance for the x-ray source ##
pitch = 0.86; % fraction of beam with couch translates per rotation

%% --- make the figure prompt for number of angles and beamlets
str = inputdlg({'Enter number of calc cores', 'Enter number of angles (51 default)', ...
    'Enter number of beamlets (64 default)'}, 'input', [1,35], {'1', '51', '64'});
num_batches = str2double(str{1}); % number of cores you want to run the beam calculation 
% -- (3 for a 4-core comp to prevent lockdown)
N_angles = str2double(str{2}); % 51 for full resolution
Mxp = str2double(str{3}); % Mxp = 64;  number of MLC leaves;
Nyp = 1;  % always 1 for Tomo due to binary mlc

% define the overall beam field size for each beam angle
% beam is 40 cm wide in transverse direction and 1-5 cm (usually 2) in y
% direction.
% isocenter is 85 cm from source, ends of jaws are 23 cm from source
xpmin = -20.0;      % -field width / 2
xpmax = 20.0;       % +field width / 2
% ypmin = -0.3125;  % total jaw width is 0.625 cm
% ypmax = 0.3125;
% ypmin = -0.5;       % total jaw width is 1 cm
% ypmax = 0.5;
ypmin = -1.25;    % total jaw width is 2.5 cm
ypmax = 1.25;
% y-prime points in the z-direction in the CT coordinate system

% ============================================= End of user-supplied inputs
% executable_path = 'C:\010-work\003_localGit\WiscPlan_v2\WiscPlanPhotonkV125\WiscPlanEXE\RyanCsphoton.x86.exe';

executable_path = mfilename('fullpath');
executable_path = [fileparts(fileparts(executable_path)), '\WiscPlanEXE\RyanCsphoton.x86.exe']

kernel_file = 'Kernels.mat';
geometry_file = fullfile(patient_dir, 'matlab_files\Geometry.mat');

load(geometry_file);
ROI_names = cellfun(@(c)c.name, Geometry.ROIS, 'UniformOutput', false);
[target_idx, okay] = listdlg('ListString', ROI_names, ...
    'SelectionMode', 'single', 'Name', 'Target Selection', ...
    'PromptString', 'Please select the target ROI for beamlet calc. ');
if okay ~= 1
    msgbox('Plan creation aborted');
    return;
end

targetMask = zeros(size(Geometry.data));
targetMask(Geometry.ROIS{target_idx}.ind) = 1;

% Grozomah - targetMask needs to get a 'double' matrix with the location of
% the target

targetMaskZ = sum(sum(targetMask,1),2);
zBow = (find(targetMaskZ>0, 1, 'first')-1)*Geometry.voxel_size(3) + Geometry.start(3) + ypmin;
zStern = (find(targetMaskZ>0, 1, 'last')+1)*Geometry.voxel_size(3) + Geometry.start(3) + ypmax;
[subi, subj, subk] = ind2sub(size(Geometry.data), Geometry.ROIS{target_idx}.ind);
iso = [Geometry.start(1)+Geometry.voxel_size(1)*mean(subi) ...
    Geometry.start(2)+Geometry.voxel_size(2)*mean(subj) 0];


% flags used to select which calculations will be set up
Condor_flag = 1;
ptvInd = target_idx;  % PTV index in Geometry.ROIS
fieldWidth = ypmax - ypmin;

% total number of rotations required for treatment
Nrot = ceil(abs(zBow - zStern)/(pitch*fieldWidth));

% Nphi = Nrot*51;  % number of angles used in the calculation
Nphi = Nrot * N_angles;  % Grozomah

% define the limits of the angles that will be used for the calculation
% ##phimin = 0;  % starting angle in radians
% ##phimax = 2*pi*Nphi;

phi = [0:Nphi-1]/Nphi *2*pi*Nrot;

%% account for beamlet shift
for scenario_i = 1:numel(OptGoals.sss_scene_list)
    patient_dir = [beamlet_dir '\scenario' num2str(scenario_i)];
    mkdir(patient_dir)
    condor_folder = patient_dir;
    winxp_folder = 'winxp';

    % create names for condor input and output folders
    input_folder = '.';
    output_folder = '.';

    % name of the convolution/superposition executable, which will be in the
    % 'code' folder of each respective run type folder
    condor_exectuable_name = 'convolutionCondor';  % relative path on the cluster where code will be
    winxp_executable_name = 'convolution.exe';
    matlab_executable_name = 'convolution_mex';  % name of the Matlab version of the dose calculation code

    % set the beam parameters, assuming a helical beam trajectory
    % folders that will be inside the 'input' folder
    beamspec_folder = 'beamspecfiles';    % directory where beam files will be stored
    beamspec_batches_folder = 'beamspecbatches';
    beamspec_batch_base_name = 'beamspecbatch';  % base name for a beamlet batch file
    kernel_folder = 'kernelfiles';  % folder where kernel information will be saved
    kernel_filenames_condor = 'kernelFilenamesCondor.txt';
    kernel_filenames_winxp = 'kernelFilenamesWinXP.txt';

    % output folders
    beamlet_batch_base_name = 'beamletbatch';  % base name for a dose batch file

    geometry_header_filename = 'geometryHeader.txt';
    geometry_density_filename = 'density.bin';  % save the density, not the Hounsfield units!

    % end of user-defined parameters

    % check the validity of the user-defined variables
    if xpmin >= xpmax
        error('xpmin must be less than xpmax.');
    end
    if ypmin >= ypmax
        error('ypmin must be less than ypmax.');b
    end

    % if phimin > phimax
    %     error('phimin must be less than or equal to phimax.');
    % end

    if Mxp <= 0 || Nyp <= 0 || Nphi <= 0
        error('Mxp, Nyp, and Nphi must be greater than zero.');
    end

    if SAD < 50
        error('It is recommended that the SAD be greater than 50 cm.');
    end

    % the xy plane is perpendicular to the isocenter axis of the linac gantry

    % size of each beam aperture, making them vectors so extension to
    % non-uniform aperture sizes becomes obvious
    del_xp = (xpmax - xpmin)/Mxp;
    del_yp = (ypmax - ypmin)/Nyp;

    % Calculate the xp and yp offsets, which lie at the centers of the
    % apertures.
    xp = [xpmin:del_xp:xpmax-del_xp] + del_xp/2;
    yp = [ypmin:del_yp:ypmax-del_yp] + del_yp/2;

    [M,N,Q] = size(Geometry.rhomw);

    START = single(Geometry.start - iso);
    INC = single(Geometry.voxel_size);

    % Grozomah ## 
    % START(1) = START(1)/10;
    % START(2) = START(2)/10;
    % INC(1) = INC(1)/10;
    % INC(2) = INC(2)/10;

    % END= START+[32,32,40].*INC

    % define the tumor mask
    tumorMask = zeros(size(Geometry.rhomw),'single');
    tumorMask(Geometry.ROIS{ptvInd}.ind) = 1;

    BW = bwdist(tumorMask);
    tumorMaskExp = tumorMask;
    tumorMaskExp(BW <= 4) = 1;

    P = zeros(Mxp,Nphi);

    fprintf('Checking beam''s eye view ...\n');
    for p=1:Nphi
        % ir and jr form the beam's eye view (BEV)
        ir = [-sin(phi(p)); cos(phi(p)); 0];
        jr = [0 0 1]';
        % kr denotes the beam direction
        kr = [cos(phi(p)); sin(phi(p)); 0];

        for m=1:Mxp
            point1 = single(-kr*SAD + [0 0 zBow + pitch*fieldWidth*phi(p)/(2*pi)]');  % source point
            point2 = single(point1 + (SAD*kr + ir*xp(m))*10);
            [indVisited,deffVisited] = singleRaytraceClean(tumorMaskExp,START,INC,point1,point2);
            if ~isempty(indVisited)
                P(m,p) = max(deffVisited);
            end
        end
    end
    fprintf('Finished checking BEV\n');

    % load data required for the dose calculator
    load(kernel_file);

    Geometry.rhomw(Geometry.rhomw < 0) = 0;
    Geometry.rhomw(Geometry.rhomw < 0.0013) = 0.0013;  % fill blank voxels with air

    % convert Geometry and kernels to single
    f = fieldnames(Kernels);
    for k=1:length(f)
        if isnumeric(getfield(Kernels,f{k}))
            Kernels = setfield(Kernels,f{k},single(getfield(Kernels,f{k})));    
        end
    end

    f = fieldnames(Geometry);
    for k=1:length(f)
        if isnumeric(getfield(Geometry,f{k}))
            Geometry = setfield(Geometry,f{k},single(getfield(Geometry,f{k})));    
        end
    end

    % account for isocenter
    Geometry.start_nominal = single(Geometry.start - iso);

    % change Condor folder names as appropriate
    
    
    
    % do the isocenter shift
    shift = OptGoals.sss_scene_list{scenario_i}; % Y X Z
    iso = [iso(1)+Geometry.voxel_size(1)*shift(1) ...
        iso(2)+Geometry.voxel_size(2)*shift(2) ...
        iso(3)+Geometry.voxel_size(3)*shift(3)];
    Geometry.start = Geometry.start_nominal- iso;

    % find the total number of beams
    Nbeam = Nphi*Mxp*Nyp;

    batch_num = 0;  % start the count for the number of total batches

    % fill up a cell array of beam structures, grouped by batch
    clear batches;
    batch_num = 0;

    batches = cell(1,Nrot);  % start the batches cell array (cell array of beam batches)
    rotNum = 0;

    % calculate beams for all source directions and apertures
    for k=1:Nphi  % loop through all gantry angles
        % calculate the source location for a helical trajectory
        beam.SAD = single(SAD);

        % the kp vector is the beam direction, ip and jp span the beam's eye view
        beam.ip = single([-sin(phi(k)) cos(phi(k)) 0]);
        beam.jp = single([0 0 1]);
        beam.kp = single([cos(phi(k)) sin(phi(k)) 0]);
        beam.y_vec = single(-beam.kp*SAD + [0 0 zBow + pitch*fieldWidth*phi(k)/(2*pi)]);

        rotNumOld = rotNum;
        rotNum = floor(k/51) + 1;  % current rotation number
        if rotNum - rotNumOld > 0
            beam_num = 0; % if the rotation number has changed, start the beam count over
        end

        for m=1:Mxp  % loop through all apertures in the xp-direction
            % calculate the beam if the tomotherapy fluence value is non-zero
            if P(m,k) > 0
                num = m + (k-1)*Mxp - 1;         % beamlet number (overall)
                beam_num = beam_num + 1;

                % set the beam aperture parameters
                beam.del_xp = single(del_xp);
                beam.del_yp = single(del_yp);
                beam.xp = single(xp(m));
                beam.yp = single(0);
                beam.num = single(num);  % record the beam number to avoid any later ambiguity

                batches{rotNum}{beam_num} = beam;
            end
        end
    end
    % merge/split batches
    all_beams = horzcat(batches{:});
    num_beams_per_batch = ceil(numel(all_beams)/num_batches);
    batches = cell(num_batches,1);
    for k = 1:(num_batches)
        beams_idx = 1+num_beams_per_batch*(k-1):num_beams_per_batch*k;
        beams_idx (beams_idx>numel(all_beams)) = [];
        batches{k} = all_beams(beams_idx);
    end
    % batches{num_batches} = all_beams(1+num_beams_per_batch*(k):end);


    % Everything else in this file is related to saving the batches in a
    % useable form.
    if Condor_flag == 1
        % delete the old submission file
        err = rmdir(fullfile(condor_folder,beamspec_batches_folder),'s');
        err = rmdir(fullfile(condor_folder,kernel_folder),'s');

        % create folders where batch information will be sent
        mkdir([condor_folder '/' input_folder '/' beamspec_batches_folder]);

        % save the kernels
        save_kernels(Kernels,[condor_folder '/' input_folder '/' kernel_folder]);
        fprintf(['Successfully saved Condor kernels to ' input_folder '/' kernel_folder '\n']);

        % create kernel filenames files
        kernel_filenames_CHTC = 'kernelFilenamesCHTC.txt';
        kernel_filenames_condor = 'kernelFilenamesCondor.txt';
        fid = fopen([condor_folder '/' input_folder '/' kernel_filenames_condor],'w');
        fid2 = fopen([condor_folder '/' input_folder '/' kernel_filenames_CHTC],'w');

        fprintf(fid,'kernel_header\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/kernel_header.txt\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'kernel_header.txt'));
        fprintf(fid2,'kernel_header\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'kernel_header.txt');

        fprintf(fid,'kernel_radii\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/radii.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'radii.bin'));
        fprintf(fid2,'kernel_radii\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'radii.bin');

        fprintf(fid,'kernel_angles\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/angles.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'angles.bin'));
        fprintf(fid2,'kernel_angles\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'angles.bin');

        fprintf(fid,'kernel_energies\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/energies.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'energies.bin'));
        fprintf(fid2,'kernel_energies\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'energies.bin');

        fprintf(fid,'kernel_primary\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/primary.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'primary.bin'));
        fprintf(fid2,'kernel_primary\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'primary.bin');  

        fprintf(fid,'kernel_first_scatter\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/first_scatter.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'first_scatter.bin'));
        fprintf(fid2,'kernel_first_scatter\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'first_scatter.bin');

        fprintf(fid,'kernel_second_scatter\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/second_scatter.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'second_scatter.bin'));
        fprintf(fid2,'kernel_second_scatter\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'second_scatter.bin');

        fprintf(fid,'kernel_multiple_scatter\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/multiple_scatter.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'multiple_scatter.bin'));
        fprintf(fid2,'kernel_multiple_scatter\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'multiple_scatter.bin');

        fprintf(fid,'kernel_brem_annih\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/brem_annih.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'brem_annih.bin'));
        fprintf(fid2,'kernel_brem_annih\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'brem_annih.bin');

        fprintf(fid,'kernel_total\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/total.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'total.bin'));
        fprintf(fid2,'kernel_total\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'total.bin');

        fprintf(fid,'kernel_fluence\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/fluence.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'fluence.bin'));
        fprintf(fid2,'kernel_fluence\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'fluence.bin');

        fprintf(fid,'kernel_mu\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/mu.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'mu.bin'));
        fprintf(fid2,'kernel_mu\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'mu.bin');

        fprintf(fid,'kernel_mu_en\n');
        % fprintf(fid,['./' input_folder '/' kernel_folder '/mu_en.bin\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,kernel_folder,'mu_en.bin'));
        fprintf(fid2,'kernel_mu_en\n');
        fprintf(fid2, '%s/%s\n', kernel_folder,'mu_en.bin');

        fclose(fid);
    end

    % name for the condor submit file that will be used
    condor_submit_file = 'convolutionSubmit.txt';

    geometry_filenames_condor = 'geometryFilenamesCondor.txt';
    geometry_filenames_CHTC = 'geometryFilenamesCHTC.txt';

    % check the geometry file to ensure that it's not in Hounsfield units
    if length(find(Geometry.rhomw > 20)) || length(find(Geometry.rhomw < 0))
        error('Double check the Geometry structure, it may still be in Hounsfield units!');
    end

    geometry_folder = 'geometryfiles';
    batch_output_folder = 'batchoutput';  % folder to which stdout will be printed
    beamlet_batches_folder = 'beamletbatches';  % folder where resulting beamlet batches will be stored

    if Condor_flag == 1
        mkdir([condor_folder '/' output_folder '/' beamlet_batches_folder]);
        mkdir([condor_folder '/' output_folder '/' batch_output_folder]);

        save_geometry(Geometry,[condor_folder '/' input_folder '/' geometry_folder],geometry_header_filename,geometry_density_filename);
        fprintf(['Successfully saved Condor geometry to ' input_folder '/' geometry_folder '\n']);

        % create geometry filenames files
        fid = fopen([condor_folder '/' input_folder '/' geometry_filenames_condor],'w');
        fid2 = fopen([condor_folder '/' input_folder '/' geometry_filenames_CHTC],'w');

        fprintf(fid,'geometry_header\n');
        % fprintf(fid,['./' input_folder '/' geometry_folder '/' geometry_header_filename '\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,geometry_folder,geometry_header_filename));
        fprintf(fid2,'geometry_header\n');
        fprintf(fid2, '%s/%s\n', geometry_folder,'geometryHeader.txt');

        fprintf(fid,'geometry_density\n');
        % fprintf(fid,['./' input_folder '/' geometry_folder '/' geometry_density_filename '\n']);
        fprintf(fid,'%s\n',fullfile(patient_dir,input_folder,geometry_folder,geometry_density_filename));
        fprintf(fid2,'geometry_density\n');
        fprintf(fid2, '%s/%s\n', geometry_folder,'density.bin');
        fclose(fid);

        % write command file
        % TODO consistent naming throughout script
        for k = 1:numel(batches)
            fid = fopen(fullfile(condor_folder,sprintf('run%d.cmd',k-1)), 'w');
            fprintf(fid, '"%s" "%s" "%s" "%s" "%s"', executable_path,...
                fullfile(patient_dir, kernel_filenames_condor),...
                fullfile(patient_dir, geometry_filenames_condor),...
                fullfile(patient_dir, 'beamspecbatches', sprintf('beamspecbatch%d.txt',k-1)),...
                fullfile(patient_dir, sprintf('batch_dose%d.bin',k-1)));
            fclose(fid);
        end

        % write the condor submit file
    %     beamspec_batch_filename = ['./' input_folder '/' beamspec_batches_folder '/' beamspec_batch_base_name '$(Process).txt'];
    %     beamlet_batch_filename = ['./' output_folder '/' beamlet_batches_folder '/' beamlet_batch_base_name '$(Process).bin'];
    %     fid = fopen([condor_folder '/' condor_submit_file],'w');
    %     fprintf(fid,'###############################################################\n');
    %     fprintf(fid,'# Condor submission script for convolution/superposition code\n');
    %     fprintf(fid,'###############################################################\n\n');
    %     fprintf(fid,'copy_to_spool = false\n');
    %     fprintf(fid,['Executable   = ' code_folder '/' condor_exectuable_name '\n']);
    %     fprintf(fid,['arguments    = ' input_folder '/' kernel_filenames_condor ' ' input_folder '/' geometry_filenames_condor ' ' beamspec_batch_filename ' ' beamlet_batch_filename '\n']);
    %     fprintf(fid,['Output       = ./' output_folder '/' batch_output_folder '/batchout$(Process).txt\n']);
    %     fprintf(fid,['Log          = ./' output_folder '/' batch_output_folder '/log.txt\n']);
    %     fprintf(fid,['Queue  ' num2str(Nrot)]);
    %     fclose(fid);



    % % write the condor submit file
    %      beamspec_batch_filename = ['./' input_folder '/' beamspec_batches_folder '/' beamspec_batch_base_name '$(Process).txt'];
    %      beamlet_batch_filename = ['./' output_folder '/' beamlet_batches_folder '/' beamlet_batch_base_name '$(Process).bin'];
         fid = fopen([condor_folder '/' condor_submit_file],'w');
         fprintf(fid,'###############################################################\n');
         fprintf(fid,'# Condor submission script for convolution/superposition code\n');
         fprintf(fid,'###############################################################\n\n');
         fprintf(fid,'copy_to_spool = false\n');
         fprintf(fid,['Executable   = ' condor_exectuable_name '\n']);
         fprintf(fid,['Arguments    = ' kernel_filenames_CHTC ' ' geometry_filenames_CHTC ' ' beamspec_batch_base_name '$(Process).txt ' 'batch_dose$(Process).bin\n']);
         fprintf(fid,['Transfer_input_files = ' kernel_folder ',' geometry_folder ',' beamspec_batches_folder '/' beamspec_batch_base_name '$(Process).txt' ',' kernel_filenames_CHTC ',' geometry_filenames_CHTC '\n']);
         fprintf(fid,['Request_memory = 1000' '\n']);
         fprintf(fid,['Request_disk = 500000' '\n']);
         fprintf(fid,['Output       = $(Cluster).out' '\n']);
         fprintf(fid,['Log          = $(Cluster).log' '\n']);
         fprintf(fid,['Error          = $(Cluster).err' '\n']);
         fprintf(fid,['Queue  ' num2str(num_batches) '\n']);
    %      fclose(fid);
    end


    % write the batches to files
    for n=1:numel(batches)
        batch = batches{n};  % current batch

        if Condor_flag == 1
            save_beamspec_batch(batch,[condor_folder '/' input_folder '/' beamspec_batches_folder],[beamspec_batch_base_name num2str(n-1) '.txt']);
        end
    end


    all_beams{1}.Mxp = Mxp;
    all_beams{1}.N_angles = N_angles;
    all_beams{1}.num_batches = num_batches;
    save([condor_folder '\all_beams.mat'], 'all_beams');
    % for k = 1:numel(batches)
    %         system([fullfile(patient_dir,sprintf('run%d.cmd',k-1)) ' &']);
    %     end 

    % Ask for User option to run the dose calculation locally on the computer 
    % or just to get necessary files for CHTC server
    % 'y' means run locally, 'n' means not to run locally on the computer

    strBeamlet = '';
    while(1)
        if strcmpi('y',strBeamlet)
            break;
        elseif strcmpi('n',strBeamlet)
            break;
        end

%         strBeamlet = input('Run beamlet batches dose calculation locally? y/n \n','s');
        strBeamlet = 'y'; %bypass question
    end

    t = datetime('now');
    disp(['Calculating ' num2str(size(all_beams, 2)) ' beamlets in ' num2str(size(batches, 1))...
        ' batches. Start: ' datestr(t)])

    if(strcmpi('y',strBeamlet)) 
        for k = 1:numel(batches)
            system([fullfile(patient_dir,sprintf('run%d.cmd',k-1)) ' &']);
        end 
    end

end % end of scenario loop

end