diff --git a/matRad_aperture2collimation.m b/matRad_aperture2collimation.m index 1ad041a53..85aef0cf7 100644 --- a/matRad_aperture2collimation.m +++ b/matRad_aperture2collimation.m @@ -51,7 +51,7 @@ pln.propStf.bixelWidth = 'field'; pln.propStf.collimation.convResolution = 0.5; %[mm] pln.propStf.collimation.fieldWidth = fieldWidth; - +pln.propStf.collimation.leafWidth = leafWidth; % %[bixelFieldX,bixelFieldY] = ndgrid(-fieldWidth/2:bixelWidth:fieldWidth/2,-fieldWidth/2:leafWidth:fieldWidth/2); diff --git a/matRad_calcParticleDoseMCtopas.m b/matRad_calcParticleDoseMCtopas.m index 57c0cc25c..ffd8f4b0d 100644 --- a/matRad_calcParticleDoseMCtopas.m +++ b/matRad_calcParticleDoseMCtopas.m @@ -24,13 +24,13 @@ % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % -% Copyright 2019 the matRad development team. -% -% This file is part of the matRad project. It is subject to the license -% terms in the LICENSE file found in the top-level directory of this -% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part -% of the matRad project, including this file, may be copied, modified, -% propagated, or distributed except according to the terms contained in the +% Copyright 2019 the matRad development team. +% +% This file is part of the matRad project. It is subject to the license +% terms in the LICENSE file found in the top-level directory of this +% distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part +% of the matRad project, including this file, may be copied, modified, +% propagated, or distributed except according to the terms contained in the % LICENSE file. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -51,48 +51,47 @@ matRad_cfg.dispWarning('Variance scoring for TOPAS not yet supported.'); end -if isfield(pln,'propMC') && isfield(pln.propMC,'config') +if isfield(pln,'propMC') && isfield(pln.propMC,'config') if isa(pln.propMC.config,'MatRad_TopasConfig') matRad_cfg.dispInfo('Using given Topas Configuration in pln.propMC.config!\n'); topasConfig = pln.propMC.config; - else + else %Create a default instance of the configuration - topasConfig = MatRad_TopasConfig(); - end - %Overwrite parameters - %mc = metaclass(topasConfig); %get metaclass information to check if we can overwrite properties - - if isstruct(pln.propMC.config) - props = fieldnames(pln.propMC.config); - for fIx = 1:numel(props) - fName = props{fIx}; - if isprop(topasConfig,fName) - %We use a try catch block to catch errors when trying - %to overwrite protected/private properties instead of a - %metaclass approach - try - topasConfig.(fName) = pln.propMC.config.(fName); - catch - matRad_cfg.dispWarning('Property ''%s'' for MatRad_TopasConfig will be omitted due to protected/private access or invalid value.',fName); - end - else - matRad_cfg.dispWarning('Unkown property ''%s'' for MatRad_TopasConfig will be omitted.',fName); + topasConfig = MatRad_TopasConfig(); + end + %Overwrite parameters + %mc = metaclass(topasConfig); %get metaclass information to check if we can overwrite properties + + if isstruct(pln.propMC.config) + props = fieldnames(pln.propMC.config); + for fIx = 1:numel(props) + fName = props{fIx}; + if isprop(topasConfig,fName) + %We use a try catch block to catch errors when trying + %to overwrite protected/private properties instead of a + %metaclass approach + try + topasConfig.(fName) = pln.propMC.config.(fName); + catch + matRad_cfg.dispWarning('Property ''%s'' for MatRad_TopasConfig will be omitted due to protected/private access or invalid value.',fName); end + else + matRad_cfg.dispWarning('Unkown property ''%s'' for MatRad_TopasConfig will be omitted.',fName); end - else - matRad_cfg.dispError('Invalid Configuration in pln.propMC.config'); end + else + matRad_cfg.dispError('Invalid Configuration in pln.propMC.config'); end else topasConfig = MatRad_TopasConfig(); end - + if ~calcDoseDirect matRad_cfg.dispError('matRad so far only supports direct dose calculation for TOPAS!\n'); end -if ~isfield(pln.propStf,'useRangeShifter') +if ~isfield(pln.propStf,'useRangeShifter') pln.propStf.useRangeShifter = false; end @@ -107,9 +106,9 @@ for s = 1:ct.numOfCtScen cubeHUresampled{s} = matRad_interp3(dij.ctGrid.x, dij.ctGrid.y', dij.ctGrid.z,ct.cubeHU{s}, ... - dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear'); + dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear'); cubeResampled{s} = matRad_interp3(dij.ctGrid.x, dij.ctGrid.y', dij.ctGrid.z,ct.cube{s}, ... - dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear'); + dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'linear'); end %Allocate temporary resampled CT @@ -150,7 +149,7 @@ % manipulate isocenter for k = 1:length(stf) stf(k).isoCenter = stf(k).isoCenter + pln.multScen.isoShift(shiftScen,:); - end + end for ctScen = 1:pln.multScen.numOfCtScen for rangeShiftScen = 1:pln.multScen.totNumRangeScen @@ -160,14 +159,14 @@ %TOPAS here) ctR.cubeHU = cubeHUresampled(ctScen); ctR.cube = cubeResampled(ctScen); - - topasConfig.writeAllFiles(ctR,pln,stf,topasBaseData,w); + + topasConfig.writeAllFiles(ctR,pln,stf,topasBaseData,w); % Run simulation for current scenario cd(topasConfig.workingDir); for beamIx = 1:numel(stf) - for runIx = 1:topasConfig.numOfRuns + for runIx = 1:topasConfig.numOfRuns fname = sprintf('%s_field%d_run%d',topasConfig.label,beamIx,runIx); topasCall = sprintf('%s %s.txt > %s.out 2> %s.err',topasConfig.topasExecCommand,fname,fname,fname); if topasConfig.parallelRuns @@ -183,7 +182,7 @@ matRad_cfg.dispError('TOPAS simulation exited with error code %d\n',status); end end - + if topasConfig.parallelRuns runsFinished = false; pause('on'); @@ -193,7 +192,7 @@ runsFinished = all(fin); end end - + end cd(currDir); @@ -201,7 +200,7 @@ %% Simulation finished - read out volume scorers from topas simulation topasCubes = matRad_readTopasData(topasConfig.workingDir); - + fnames = fieldnames(topasCubes); dij.MC_tallies = fnames; for f = 1:numel(fnames) diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m index fc56ffbdf..ee05e39bc 100644 --- a/matRad_calcPhotonDoseMCtopas.m +++ b/matRad_calcPhotonDoseMCtopas.m @@ -126,14 +126,18 @@ topasConfig.numOfRuns = 1; %matRad_cfg.propMC.topas_defaultNumBatches; topasConfig.beamProfile = "uniform"; %"virtualGaussian"; %Collect weights -w = zeros(sum([stf(:).totalNumOfBixels]),1); -ct = 1; -for i = 1:length(stf) - for j = 1:stf(i).numOfRays - rayBix = stf(i).numOfBixelsPerRay(j); - w(ct:ct+rayBix-1) = stf(i).ray(j).weight; - ct = ct + rayBix; +if calcDoseDirect + w = zeros(sum([stf(:).totalNumOfBixels]),1); + ct = 1; + for i = 1:length(stf) + for j = 1:stf(i).numOfRays + rayBix = stf(i).numOfBixelsPerRay(j); + w(ct:ct+rayBix-1) = stf(i).ray(j).weight; + ct = ct + rayBix; + end end +else + w = ones(sum([stf(:).totalNumOfBixels]),1); end currDir = cd; diff --git a/topas/MatRad_TopasConfig.m b/topas/MatRad_TopasConfig.m index a8f074c10..56f3514ef 100644 --- a/topas/MatRad_TopasConfig.m +++ b/topas/MatRad_TopasConfig.m @@ -110,8 +110,7 @@ %Let's set some default commands taken from topas installation %instructions for mac & debain/ubuntu if ispc %We assume topas is installed in wsl (since no windows version) - %obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; ~/topas/bin/topas'; - obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; export LD_LIBRARY_PATH=~/topas/libexternal/:$LD_LIBRARY_PATH; ~/topas/topas'; %<-- for Pia's laptop uncomment + obj.topasExecCommand = 'wsl export TOPAS_G4_DATA_DIR=~/G4Data; ~/topas/bin/topas'; elseif ismac obj.topasExecCommand = 'export TOPAS_G4_DATA_DIR=/Applications/G4Data; export QT_QPA_PLATFORM_PLUGIN_PATH=/Applications/topas/Frameworks; /Applications/topas/bin/topas'; elseif isunix @@ -673,7 +672,7 @@ function writePatient(obj,ct,pln) end %cubeExport = obj.materialConversion; %'RSP'; %'HUSchneiderToWater'; 'CustomWaterRSP' - cubeExport = 'CustomWaterRSP'; + cubeExport = 'RSP'; rspCubeMethod = 2; checkMaterial = false; diff --git a/topas/TOPAS_beamSetup_mlc.txt.in b/topas/TOPAS_beamSetup_mlc.txt.in index 0ba1aed36..4cdfa44fa 100644 --- a/topas/TOPAS_beamSetup_mlc.txt.in +++ b/topas/TOPAS_beamSetup_mlc.txt.in @@ -3,7 +3,7 @@ s:Ge/MultiLeafCollimatorA/Parent = "Nozzle" s:Ge/MultiLeafCollimatorA/Material = "Aluminum" d:Ge/MultiLeafCollimatorA/TransX = 0.0 cm d:Ge/MultiLeafCollimatorA/TransY = 0.0 cm -d:Ge/MultiLeafCollimatorA/TransZ = 0.0 cm +d:Ge/MultiLeafCollimatorA/TransZ = 7.5 cm d:Ge/MultiLeafCollimatorA/RotX = 0.0 deg d:Ge/MultiLeafCollimatorA/RotY = 0.0 deg d:Ge/MultiLeafCollimatorA/RotZ = 90.0 deg diff --git a/topas/matRad_TopasConfigPhotons.m b/topas/matRad_TopasConfigPhotons.m index d4d16180c..ec9a3ce06 100644 --- a/topas/matRad_TopasConfigPhotons.m +++ b/topas/matRad_TopasConfigPhotons.m @@ -430,42 +430,42 @@ function writeStfFields(obj,ct,stf,pln,baseData,w) if exist('TOPAS_mlcSetup') fprintf(fileID,'%s\n',TOPAS_mlcSetup); - %Erstmal nur für einen ray! - [l,t]=size([stf(1).ray.shapes(:).leftLeafPos]); %there are l leaves and t times/shapes + %Erstmal nur f??r einen ray! + [numOfLeaves,leafTimes]=size([stf(1).ray.shapes(:).leftLeafPos]); %there are #numOfLeaves leaves and #leafTimes times/shapes leftLeafPos = [stf(1).ray.shapes(:).leftLeafPos]; rightLeafPos = [stf(1).ray.shapes(:).rightLeafPos]; %Set MLC paramters as in TOPAS example file https://topas.readthedocs.io/en/latest/parameters/geometry/specialized.html#multi-leaf-collimator fprintf(fileID,'d:Ge/MultiLeafCollimatorA/MaximumLeafOpen = %f cm\n',15); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Thickness = %f cm\n',15); fprintf(fileID,'d:Ge/MultiLeafCollimatorA/Length = %f cm\n',6); - fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', l) - fprintf(fileID,num2str(ones(1,l),' % 2d')); - fprintf(fileID,' cm \n'); - fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',l); - for i=1:l + fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/Widths = %i ', numOfLeaves); + fprintf(fileID,num2str(pln.propStf.collimation.leafWidth*ones(1,numOfLeaves),' % 2d')); + fprintf(fileID,' mm \n'); + fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XPlusLeavesOpen = %i ',numOfLeaves); + for i=1:numOfLeaves fprintf( fileID,'Tf/LeafXPlus%i/Value ',i); end fprintf(fileID,'mm \n'); - fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',l); - for i=1:l + fprintf(fileID,'dv:Ge/MultiLeafCollimatorA/XMinusLeavesOpen = %i ',numOfLeaves); + for i=1:numOfLeaves fprintf( fileID,'Tf/LeafXMinus%i/Value ',i); end fprintf(fileID,'mm \n'); - for i=1:l + for i=1:numOfLeaves fprintf(fileID,'s:Tf/LeafXMinus%i/Function = "Step"\n',i); - fprintf(fileID,'dv:Tf/LeafXMinus%i/Times = %i ', i,t); - fprintf(fileID,num2str([1:t]*10,' % 2d')); + fprintf(fileID,'dv:Tf/LeafXMinus%i/Times = %i ', i,leafTimes); + fprintf(fileID,num2str([1:leafTimes]*10,' % 2d')); fprintf(fileID,' ms\n'); - fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,t); + fprintf(fileID,'dv:Tf/LeafXMinus%i/Values = %i ', i,leafTimes); fprintf(fileID,num2str(leftLeafPos(i,:),' % 2d')); fprintf(fileID,' mm\n\n'); fprintf(fileID,'s:Tf/LeafXPlus%i/Function = "Step"\n',i); - fprintf(fileID,'dv:Tf/LeafXPlus%i/Times = %i ',i,t); - fprintf(fileID,num2str([1:t]*10,' % 2d')); + fprintf(fileID,'dv:Tf/LeafXPlus%i/Times = %i ',i,leafTimes); + fprintf(fileID,num2str([1:leafTimes]*10,' % 2d')); fprintf(fileID,' ms\n'); - fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,t); + fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes); fprintf(fileID,num2str(rightLeafPos(i,:),' % 2d')); fprintf(fileID,' mm\n\n');