diff --git a/MatRad_Config.m b/MatRad_Config.m index 974432f1e..c4403a72e 100644 --- a/MatRad_Config.m +++ b/MatRad_Config.m @@ -186,7 +186,7 @@ function setDefaultProperties(obj) % Set default histories for MonteCarlo here if necessary % obj.propMC.defaultNumHistories = 100; - %obj.propMC.default_photon_engine = 'ompMC'; + obj.propMC.default_photon_engine = 'ompMC'; obj.propMC.default_proton_engine = 'MatRad_MCsquareConfig'; obj.propMC.default_carbon_engine = 'MatRad_TopasConfig'; @@ -239,7 +239,7 @@ function setDefaultPropertiesForTesting(obj) % Set default histories for MonteCarlo obj.propMC.defaultNumHistories = 100; - %obj.propMC.default_photon_engine = 'ompMC'; + obj.propMC.default_photon_engine = 'ompMC'; obj.propMC.default_proton_engine = 'MatRad_MCsquareConfig'; obj.propMC.default_carbon_engine = 'MatRad_TopasConfig'; @@ -402,6 +402,8 @@ function getEnvironment(obj) else if isfield(pln,'radiationMode') && ~isempty(pln.radiationMode) switch pln.radiationMode + case 'photons' + configName = obj.propMC.default_photon_engine; case 'protons' configName = obj.propMC.default_proton_engine; otherwise diff --git a/matRad_aperture2collimation.m b/matRad_aperture2collimation.m new file mode 100644 index 000000000..2cb30a6c7 --- /dev/null +++ b/matRad_aperture2collimation.m @@ -0,0 +1,130 @@ +function [pln,stf] = matRad_aperture2collimation(pln,stf,sequencing,apertureInfo) +% matRad function to convert sequencing information / aperture information +% into collimation information in pln and stf for field-based dose +% calculation +% +% call +% [pln,stf] = matRad_aperture2collimation(pln,stf,sequencing,apertureInfo) +% +% input +% pln: pln file used to generate the sequenced plan +% stf: stf file used to generate the sequenced plan +% sequencing: sequencing information (from resultGUI) +% apertureInfo: apertureInfo (from resultGUI) +% +% output +% pln: matRad pln struct with collimation information +% stf: matRad stf struct with shapes instead of beamlets +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2022 the matRad development team. +% Author: wahln +% +% 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%Debug visualization +visBool = false; + +bixelWidth = apertureInfo.bixelWidth; +leafWidth = bixelWidth; +convResolution = 0.5; %[mm] + +%The collimator limits are infered here from the apertureInfo. This could +%be handled differently by explicitly storing collimator info in the base +%data? +symmetricMLClimits = vertcat(apertureInfo.beam.MLCWindow); +symmetricMLClimits = max(abs(symmetricMLClimits)); +fieldWidth = 2*max(symmetricMLClimits); + +%modify basic pln variables +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); +[convFieldX,convFieldY] = meshgrid(-fieldWidth/2:convResolution:fieldWidth/2); + +%TODO: Not used in calcPhotonDose but imported from DICOM +%pln.propStf.collimation.Devices ... +%pln.propStf.collimation.numOfFields +%pln.propStf.collimation.beamMeterset + +for iBeam = 1:numel(stf) + stfTmp = stf(iBeam); + beamSequencing = sequencing.beam(iBeam); + beamAperture = apertureInfo.beam(iBeam); + + stfTmp.bixelWidth = 'field'; + + nShapes = beamSequencing.numOfShapes; + + stfTmp.numOfRays = 1;% + stfTmp.numOfBixelsPerRay = nShapes; + stfTmp.totalNumOfBixels = nShapes; + + ray = struct(); + ray.rayPos_bev = [0 0 0]; + ray.targetPoint_bev = [0 stfTmp.SAD 0]; + ray.weight = 1; + ray.energy = stfTmp.ray(1).energy; + ray.beamletCornersAtIso = stfTmp.ray(1).beamletCornersAtIso; + ray.rayCorners_SCD = stfTmp.ray(1).rayCorners_SCD; + + %ray.shape = beamSequencing.sum; + shapeTotalF = zeros(size(convFieldX)); + + ray.shapes = struct(); + for iShape = 1:nShapes + currShape = beamAperture.shape(iShape); + activeLeafPairPosY = beamAperture.leafPairPos; + F = zeros(size(convFieldX)); + if visBool + hF = figure; imagesc(F); title(sprintf('Beam %d, Shape %d',iBeam,iShape)); hold on; + end + for iLeafPair = 1:numel(activeLeafPairPosY) + posY = activeLeafPairPosY(iLeafPair); + ixY = convFieldY >= posY-leafWidth/2 & convFieldY < posY + leafWidth/2; + ixX = convFieldX >= currShape.leftLeafPos(iLeafPair) & convFieldX < currShape.rightLeafPos(iLeafPair); + ix = ixX & ixY; + F(ix) = 1; + if visBool + figure(hF); imagesc(F); drawnow; pause(0.1); + end + end + + if visBool + pause(1); close(hF); + end + + F = F*currShape.weight; + shapeTotalF = shapeTotalF + F; + + ray.shapes(iShape).convFluence = F; + ray.shapes(iShape).shapeMap = currShape.shapeMap; + ray.shapes(iShape).weight = currShape.weight; + ray.shapes(iShape).leftLeafPos = currShape.leftLeafPos; + ray.shapes(iShape).rightLeafPos = currShape.rightLeafPos; + ray.shapes(iShape).leafPairCenterPos = activeLeafPairPosY; + end + + ray.shape = shapeTotalF; + ray.weight = ones(1,nShapes); + stfTmp.ray = ray; + + stf(iBeam) = stfTmp; +end + + diff --git a/matRad_calcPhotonDoseMC.m b/matRad_calcPhotonDoseMC.m index 67e1cb721..6e21f097c 100644 --- a/matRad_calcPhotonDoseMC.m +++ b/matRad_calcPhotonDoseMC.m @@ -1,363 +1,59 @@ -function dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,nCasePerBixel,visBool) -% matRad ompMC monte carlo photon dose calculation wrapper +function dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,calcDoseDirect) +% matRad Monte Carlo photon dose calculation wrapper +% Will call the appropriate subfunction for the respective +% MC dose-calculation engine % % call -% dij = matRad_calcPhotonDoseMc(ct,stf,pln,cst,visBool) +% dij = matRad_calcParticleDoseMc(ct,stf,pln,cst) +% dij = matRad_calcParticleDoseMc(ct,stf,pln,cst,calcDoseDirect) % % input % ct: matRad ct struct % stf: matRad steering information struct % pln: matRad plan meta information struct % cst: matRad cst struct -% visBool: binary switch to enable visualization +% nHistories: number of histories per beamlet +% calcDoseDirect: (optional) binary switch to enable forward +% dose calcualtion (default: false) % output % dij: matRad dij struct % % References -% - % -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% - % -% Copyright 2018 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 2020 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. % -% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -matRad_cfg = MatRad_Config.instance(); +matRad_cfg = MatRad_Config.instance(); -% disable visualiazation by default if nargin < 6 - visBool = false; -end - -if nargin < 5 - nCasePerBixel = matRad_cfg.propMC.ompMC_defaultHistories; - matRad_cfg.dispInfo('Using default number of Histories per Bixel: %d\n',nCasePerBixel); -end - -fileFolder = fileparts(mfilename('fullpath')); - -if ~matRad_checkMexFileExists('omc_matrad') %exist('matRad_ompInterface','file') ~= 3 - matRad_cfg.dispWarning('Compiled mex interface not found. Trying to compile the ompMC interface on the fly!'); - try - matRad_compileOmpMCInterface(); - catch MException - matRad_cfg.dispError('Could not find/generate mex interface for MC dose calculation.\nCause of error:\n%s\n Please compile it yourself (preferably with OpenMP support).',MException.message); - end -end - -matRad_calcDoseInit; - -% gaussian filter to model penumbra from (measured) machine output / see diploma thesis siggel 4.1.2 -if isfield(machine.data,'penumbraFWHMatIso') - penumbraFWHM = machine.data.penumbraFWHMatIso; -else - penumbraFWHM = 5; - matRad_cfg.dispWarning('photon machine file does not contain measured penumbra width in machine.data.penumbraFWHMatIso. Assuming 5 mm.'); -end - -sourceFWHM = penumbraFWHM * machine.meta.SCD/(machine.meta.SAD - machine.meta.SCD); -sigmaGauss = sourceFWHM / sqrt(8*log(2)); % [mm] - -% set up arrays for book keeping -dij.bixelNum = NaN*ones(dij.totalNumOfBixels,1); -dij.rayNum = NaN*ones(dij.totalNumOfBixels,1); -dij.beamNum = NaN*ones(dij.totalNumOfBixels,1); - -dij.numHistoriesPerBeamlet = nCasePerBixel; - -omcFolder = [matRad_cfg.matRadRoot filesep 'ompMC']; -%omcFolder = [matRad_cfg.matRadRoot filesep 'submodules' filesep 'ompMC']; - -%% Setup OmpMC options / parameters - -%display options -ompMCoptions.verbose = matRad_cfg.logLevel - 1; - -% start MC control -ompMCoptions.nHistories = nCasePerBixel; -ompMCoptions.nSplit = 20; -ompMCoptions.nBatches = 10; -ompMCoptions.randomSeeds = [97 33]; - -%start source definition -ompMCoptions.spectrumFile = [omcFolder filesep 'spectra' filesep 'mohan6.spectrum']; -ompMCoptions.monoEnergy = 0.1; -ompMCoptions.charge = 0; -ompMCoptions.sourceGeometry = 'gaussian'; -ompMCoptions.sourceGaussianWidth = 0.1*sigmaGauss; - -% start MC transport -ompMCoptions.dataFolder = [omcFolder filesep 'data' filesep]; -ompMCoptions.pegsFile = [omcFolder filesep 'pegs4' filesep '700icru.pegs4dat']; -ompMCoptions.pgs4formFile = [omcFolder filesep 'pegs4' filesep 'pgs4form.dat']; - -ompMCoptions.global_ecut = 0.7; -ompMCoptions.global_pcut = 0.010; - -% Relative Threshold for dose -ompMCoptions.relDoseThreshold = 1 - matRad_cfg.propDoseCalc.defaultLateralCutOff; - -% Output folders -ompMCoptions.outputFolder = [omcFolder filesep 'output' filesep]; - -% Create Material Density Cube -material = cell(4,5); -material{1,1} = 'AIR700ICRU'; -material{1,2} = -1024; -material{1,3} = -974; -material{1,4} = 0.001; -material{1,5} = 0.044; -material{2,1} = 'LUNG700ICRU'; -material{2,2} = -974; -material{2,3} = -724; -material{2,4} = 0.044; -material{2,5} = 0.302; -material{3,1} = 'ICRUTISSUE700ICRU'; -material{3,2} = -724; -material{3,3} = 101; -material{3,4} = 0.302; -material{3,5} = 1.101; -material{4,1} = 'ICRPBONE700ICRU'; -material{4,2} = 101; -material{4,3} = 1976; -material{4,4} = 1.101; -material{4,5} = 2.088; - -% conversion from HU to densities & materials -for s = 1:ct.numOfCtScen - - HUcube{s} = matRad_interp3(dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z,ct.cubeHU{s}, ... - dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'); - - % projecting out of bounds HU values where necessary - if max(HUcube{s}(:)) > material{end,3} - matRad_cfg.dispWarning('Projecting out of range HU values'); - HUcube{s}(HUcube{s}(:) > material{end,3}) = material{end,3}; - end - if min(HUcube{s}(:)) < material{1,2} - matRad_cfg.dispWarning('Projecting out of range HU values'); - HUcube{s}(HUcube{s}(:) < material{1,2}) = material{1,2}; - end - - % find material index - cubeMatIx{s} = NaN*ones(dij.doseGrid.dimensions,'int32'); - for i = size(material,1):-1:1 - cubeMatIx{s}(HUcube{s} <= material{i,3}) = i; - end - - % create an artificial HU lookup table - hlut = []; - for i = 1:size(material,1) - hlut = [hlut;material{i,2} material{i,4};material{i,3}-1e-10 material{i,5}]; % add eps for interpolation - end - - cubeRho{s} = interp1(hlut(:,1),hlut(:,2),HUcube{s}); - + calcDoseDirect = false; end -ompMCgeo.material = material; +% load appropriate config from pln or from class +pln = matRad_cfg.getDefaultClass(pln,'propMC'); -scale = 10; % to convert to cm +% load default parameters in case they haven't been set yet +pln = matRad_cfg.getDefaultProperties(pln,'propDoseCalc'); -ompMCgeo.xBounds = (dij.doseGrid.resolution.y * (0.5 + [0:dij.doseGrid.dimensions(1)])) ./ scale; -ompMCgeo.yBounds = (dij.doseGrid.resolution.x * (0.5 + [0:dij.doseGrid.dimensions(2)])) ./ scale; -ompMCgeo.zBounds = (dij.doseGrid.resolution.z * (0.5 + [0:dij.doseGrid.dimensions(3)])) ./ scale; - -%% debug visualization -if visBool - - figure - hold on - - axis equal - - % ct box - ctCorner1 = [ompMCgeo.xBounds(1) ompMCgeo.yBounds(1) ompMCgeo.zBounds(1)]; - ctCorner2 = [ompMCgeo.xBounds(end) ompMCgeo.yBounds(end) ompMCgeo.zBounds(end)]; - plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner1(3)],'k' ) - plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) - plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) - plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) - plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner2(3) ctCorner2(3)],'k' ) - plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) - plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) - plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) - plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' ) - plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' ) - plot3([ctCorner1(1) ctCorner1(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' ) - plot3([ctCorner2(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' ) - - xlabel('x [cm]') - ylabel('y [cm]') - zlabel('z [cm]') - - rotate3d on - +switch pln.propMC.engine + case 'ompMC' + dij = matRad_calcPhotonDoseOmpMC(ct,stf,pln,cst,nHistories,calcDoseDirect); + case 'TOPAS' + dij = matRad_calcPhotonDoseMCtopas(ct,stf,pln,cst,calcDoseDirect); + otherwise + matRad_cfg.dispError('MC engine %s not known/supported!',engine); end -%ompMC for matRad returns dose/history * nHistories. -% This factor calibrates to 1 Gy in a %(5x5)cm^2 open field (1 bixel) at -% 5cm depth for SSD = 900 which corresponds to the calibration for the -% analytical base data. -absCalibrationFactor = 3.49056 * 1e12; %Approximate! - -%Now we have to calibrate to the the beamlet width. -absCalibrationFactor = absCalibrationFactor * (pln.propStf.bixelWidth/50)^2; - -outputVariance = matRad_cfg.propMC.ompMC_defaultOutputVariance; - -if isfield(pln,'propMC') && isfield(pln.propMC,'outputVariance') - outputVariance = pln.propMC.outputVariance; -end - -%% Create beamlet source -useCornersSCD = true; %false -> use ISO corners - -scenCount = 0; - -for shiftScen = 1:pln.multScen.totNumShiftScen - - % manipulate isocenter - for k = 1:length(stf) - stf(k).isoCenter = stf(k).isoCenter + pln.multScen.isoShift(shiftScen,:); - end - - for ctScen = 1:pln.multScen.numOfCtScen - for rangeShiftScen = 1:pln.multScen.totNumRangeScen - if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) - scenCount = scenCount + 1; - - numOfBixels = [stf(:).numOfRays]; - beamSource = zeros(dij.numOfBeams, 3); - - bixelCorner = zeros(dij.totalNumOfBixels,3); - bixelSide1 = zeros(dij.totalNumOfBixels,3); - bixelSide2 = zeros(dij.totalNumOfBixels,3); - - counter = 0; - - for i = 1:dij.numOfBeams % loop over all beams - - % define beam source in physical coordinate system in cm - beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10; - - for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray! - - counter = counter + 1; - - dij.beamNum(counter) = i; - dij.rayNum(counter) = j; - dij.bixelNum(counter) = j; - - if useCornersSCD - beamletCorners = stf(i).ray(j).rayCorners_SCD; - else - beamletCorners = stf(i).ray(j).beamletCornersAtIso; - end - - % get bixel corner and delimiting vectors. - % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm - currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ scale; - bixelCorner(counter,:) = currCorner; - bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ scale - currCorner; - bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ scale - currCorner; - - if visBool - for k = 1:4 - currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10; - % rays connecting source and ray corner - plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b') - % connection between corners - lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10; - plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r') - end - end - - end - - end - - ompMCsource.nBeams = dij.numOfBeams; - ompMCsource.iBeam = dij.beamNum(:); - - % Switch x and y directions to match ompMC cs. - ompMCsource.xSource = beamSource(:,2); - ompMCsource.ySource = beamSource(:,1); - ompMCsource.zSource = beamSource(:,3); - - ompMCsource.nBixels = sum(numOfBixels(:)); - ompMCsource.xCorner = bixelCorner(:,2); - ompMCsource.yCorner = bixelCorner(:,1); - ompMCsource.zCorner = bixelCorner(:,3); - - ompMCsource.xSide1 = bixelSide1(:,2); - ompMCsource.ySide1 = bixelSide1(:,1); - ompMCsource.zSide1 = bixelSide1(:,3); - - ompMCsource.xSide2 = bixelSide2(:,2); - ompMCsource.ySide2 = bixelSide2(:,1); - ompMCsource.zSide2 = bixelSide2(:,3); - - if visBool - plot3(ompMCsource.ySource,ompMCsource.xSource,ompMCsource.zSource,'rx') - end - - %% Call the OmpMC interface - if pln.multScen.totNumScen > 1 - matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation... \n'); - else - matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation for scenario %d of %d... \n',scenCount,pln.multScen.totNumScen); - end - - ompMCgeo.isoCenter = [stf(:).isoCenter]; - - %Call the Monte Carlo simulation and catch possible mex - %interface issues - try - %If we ask for variance, a field in the dij will be filled - if outputVariance - [dij.physicalDose{ctScen,shiftScen,rangeShiftScen},dij.physicalDose_MCvar{ctScen,shiftScen,rangeShiftScen}] = omc_matrad(cubeRho{ctScen},cubeMatIx{ctScen},ompMCgeo,ompMCsource,ompMCoptions); - else - [dij.physicalDose{ctScen,shiftScen,rangeShiftScen}] = omc_matrad(cubeRho{ctScen},cubeMatIx{ctScen},ompMCgeo,ompMCsource,ompMCoptions); - end - catch ME - errorString = [ME.message '\nThis error was thrown by the MEX-interface of ompMC.\nMex interfaces can raise compatability issues which may be resolved by compiling them by hand directly on your particular system.']; - matRad_cfg.dispError(errorString); - end - - dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} * absCalibrationFactor; - if isfield(dij,'physicalDose_MCvar') - dij.physicalDose_MCvar{s} = dij.physicalDose_MCvar{ctScen,shiftScen,rangeShiftScen} * absCalibrationFactor^2; - end - - ompMCgeo.isoCenter = [stf(:).isoCenter]; - - matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n'); - - try - % wait 0.1s for closing all waitbars - allWaitBarFigures = findall(0,'type','figure','tag','TMWWaitbar'); - delete(allWaitBarFigures); - pause(0.1); - catch - end - - % manipulate isocenter back - for k = 1:length(stf) - stf(k).isoCenter = stf(k).isoCenter - pln.multScen.isoShift(shiftScen,:); - end - end - end - end -end - - end diff --git a/matRad_calcPhotonDoseMCtopas.m b/matRad_calcPhotonDoseMCtopas.m new file mode 100644 index 000000000..25331b09b --- /dev/null +++ b/matRad_calcPhotonDoseMCtopas.m @@ -0,0 +1,164 @@ +function dij = matRad_calcPhotonDoseMCtopas(ct,stf,pln,cst,calcDoseDirect) +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% matRad TOPAS Monte Carlo proton dose calculation wrapper +% This calls a TOPAS installation (not included in matRad due to +% licensing model of TOPAS) for MC simulation +% +% call +% dij = matRad_calcParticleDoseMCtopas(ct,stf,pln,cst,calcDoseDirect) +% +% input +% ct: matRad ct struct +% stf: matRad steering information struct +% pln: matRad plan meta information struct +% cst: matRad cst struct +% calcDoseDirect: binary switch to enable forward dose +% calcualtion +% output +% dij: matRad dij struct +% +% References +% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +matRad_cfg = MatRad_Config.instance(); + +if nargin < 6 + calcDoseDirect = false; +end + +% Set parameters for full Dij calculation +if ~calcDoseDirect + pln.propMC.scorer.calcDij = true; + pln.propMC.numOfRuns = 1; + + % Load class variables in pln + % for calcDoseDirect, this is already done in superior function + pln = matRad_cfg.getDefaultClass(pln,'propMC','MatRad_TopasConfig'); +else + if ~isa(pln.propMC,'MatRad_TopasConfig') + matRad_cfg.dispError('Run calcParticleDoseMCtopas through calcDoseDirectMC'); + end +end + +% load default parameters for doseCalc in case they haven't been set yet +pln = matRad_cfg.getDefaultProperties(pln,'propDoseCalc'); + +% set nested folder structure if external calculation is turned on (this will put new simulations in subfolders) +if pln.propMC.externalCalculation + pln.propMC.workingDir = [pln.propMC.thisFolder filesep 'MCrun' filesep]; + pln.propMC.workingDir = [pln.propMC.workingDir pln.radiationMode,'_',pln.machine,'_',datestr(now, 'dd-mm-yy')]; +end + +%% Initialize dose Grid as usual +matRad_calcDoseInit; + +% for TOPAS we explicitly downsample the ct to the dose grid (might not be necessary in future versions with separated grids) +[ctR,~,~] = matRad_resampleCTtoGrid(ct,cst,pln,stf); + +% overwrite CT grid in dij in case of modulation. +if isfield(ctR,'ctGrid') + dij.ctGrid = ctR.ctGrid; +end + +%% sending data to topas + +load([pln.radiationMode,'_',pln.machine]); + +pln.propMC.numOfRuns = 1; %matRad_cfg.propMC.topas_defaultNumBatches; +pln.propMC.beamProfile = 'phasespace'; %'uniform'; 'virtualGaussian'; +%Collect weights +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; + +for shiftScen = 1:pln.multScen.totNumShiftScen + + % manipulate isocenter + for k = 1:length(stf) + stf(k).isoCenter = stf(k).isoCenter + pln.multScen.isoShift(shiftScen,:); + end + + for ctScen = 1:pln.multScen.numOfCtScen + for rangeShiftScen = 1:pln.multScen.totNumRangeScen + if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) + + pln.propMC.writeAllFiles(ctR,cst,pln,stf,machine,w); + + % Run simulation for current scenario + cd(pln.propMC.workingDir); + + for beamIx = 1:numel(stf) + for runIx = 1:pln.propMC.numOfRuns + fname = sprintf('%s_field%d_run%d',pln.propMC.label,beamIx,runIx); + topasCall = sprintf('%s %s.txt > %s.out 2> %s.err',pln.propMC.topasExecCommand,fname,fname,fname); + if pln.propMC.parallelRuns + finishedFiles{runIx} = sprintf('%s.finished',fname); + delete(finishedFiles{runIx}); + topasCall = [topasCall '; touch ' finishedFiles{runIx} ' &']; + end + matRad_cfg.dispInfo('Calling TOPAS: %s\n',topasCall); + [status,cmdout] = system(topasCall,'-echo'); + if status == 0 + matRad_cfg.dispInfo('TOPAS simulation completed succesfully\n'); + else + matRad_cfg.dispError('TOPAS simulation exited with error code %d\n',status); + end + end + + if pln.propMC.parallelRuns + runsFinished = false; + pause('on'); + while ~runsFinished + pause(1); + fin = cellfun(@(f) exist(f,'file'),finishedFiles); + runsFinished = all(fin); + end + end + + end + + cd(currDir); + + %% Simulation finished - read out volume scorers from topas simulation + + topasCubes = matRad_readTopasData(pln.propMC.workingDir); + + fnames = fieldnames(topasCubes); + dij.MC_tallies = fnames; + for f = 1:numel(fnames) + dij.(fnames{f}){1} = sum(w)*reshape(topasCubes.(fnames{f}),[],1); + end + end + end + end + + % manipulate isocenter back + for k = 1:length(stf) + stf(k).isoCenter = stf(k).isoCenter - pln.multScen.isoShift(shiftScen,:); + end +end diff --git a/matRad_calcPhotonDoseOmpMC.m b/matRad_calcPhotonDoseOmpMC.m new file mode 100644 index 000000000..96c83e88f --- /dev/null +++ b/matRad_calcPhotonDoseOmpMC.m @@ -0,0 +1,371 @@ +function dij = matRad_calcPhotonDoseMC(ct,stf,pln,cst,nCasePerBixel,calcDoseDirect,visBool) +% matRad ompMC monte carlo photon dose calculation wrapper +% +% call +% dij = matRad_calcPhotonDoseMc(ct,stf,pln,cst,visBool) +% +% input +% ct: matRad ct struct +% stf: matRad steering information struct +% pln: matRad plan meta information struct +% cst: matRad cst struct +% visBool: binary switch to enable visualization +% output +% dij: matRad dij struct +% +% References +% - +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Copyright 2018 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. +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +matRad_cfg = MatRad_Config.instance(); + +if nargin < 6 + calcDoseDirect = false; +end + +if calcDoseDirect + matRad_cfg.dispWarning('ompMC does not support forward dose calculation natively! Will compute full dij and then apply fluence vector.'); +end + +% disable visualiazation by default +if nargin < 7 + visBool = false; +end + +if nargin < 5 + nCasePerBixel = matRad_cfg.propMC.ompMC_defaultHistories; + matRad_cfg.dispInfo('Using default number of Histories per Bixel: %d\n',nCasePerBixel); +end + +fileFolder = fileparts(mfilename('fullpath')); + +if ~matRad_checkMexFileExists('omc_matrad') %exist('matRad_ompInterface','file') ~= 3 + matRad_cfg.dispWarning('Compiled mex interface not found. Trying to compile the ompMC interface on the fly!'); + try + matRad_compileOmpMCInterface(); + catch MException + matRad_cfg.dispError('Could not find/generate mex interface for MC dose calculation.\nCause of error:\n%s\n Please compile it yourself (preferably with OpenMP support).',MException.message); + end +end + +matRad_calcDoseInit; + +% gaussian filter to model penumbra from (measured) machine output / see diploma thesis siggel 4.1.2 +if isfield(machine.data,'penumbraFWHMatIso') + penumbraFWHM = machine.data.penumbraFWHMatIso; +else + penumbraFWHM = 5; + matRad_cfg.dispWarning('photon machine file does not contain measured penumbra width in machine.data.penumbraFWHMatIso. Assuming 5 mm.'); +end + +sourceFWHM = penumbraFWHM * machine.meta.SCD/(machine.meta.SAD - machine.meta.SCD); +sigmaGauss = sourceFWHM / sqrt(8*log(2)); % [mm] + +% set up arrays for book keeping +dij.bixelNum = NaN*ones(dij.totalNumOfBixels,1); +dij.rayNum = NaN*ones(dij.totalNumOfBixels,1); +dij.beamNum = NaN*ones(dij.totalNumOfBixels,1); + +dij.numHistoriesPerBeamlet = nCasePerBixel; + +omcFolder = [matRad_cfg.matRadRoot filesep 'ompMC']; +%omcFolder = [matRad_cfg.matRadRoot filesep 'submodules' filesep 'ompMC']; + +%% Setup OmpMC options / parameters + +%display options +ompMCoptions.verbose = matRad_cfg.logLevel - 1; + +% start MC control +ompMCoptions.nHistories = nCasePerBixel; +ompMCoptions.nSplit = 20; +ompMCoptions.nBatches = 10; +ompMCoptions.randomSeeds = [97 33]; + +%start source definition +ompMCoptions.spectrumFile = [omcFolder filesep 'spectra' filesep 'mohan6.spectrum']; +ompMCoptions.monoEnergy = 0.1; +ompMCoptions.charge = 0; +ompMCoptions.sourceGeometry = 'gaussian'; +ompMCoptions.sourceGaussianWidth = 0.1*sigmaGauss; + +% start MC transport +ompMCoptions.dataFolder = [omcFolder filesep 'data' filesep]; +ompMCoptions.pegsFile = [omcFolder filesep 'pegs4' filesep '700icru.pegs4dat']; +ompMCoptions.pgs4formFile = [omcFolder filesep 'pegs4' filesep 'pgs4form.dat']; + +ompMCoptions.global_ecut = 0.7; +ompMCoptions.global_pcut = 0.010; + +% Relative Threshold for dose +ompMCoptions.relDoseThreshold = 1 - matRad_cfg.propDoseCalc.defaultLateralCutOff; + +% Output folders +ompMCoptions.outputFolder = [omcFolder filesep 'output' filesep]; + +% Create Material Density Cube +material = cell(4,5); +material{1,1} = 'AIR700ICRU'; +material{1,2} = -1024; +material{1,3} = -974; +material{1,4} = 0.001; +material{1,5} = 0.044; +material{2,1} = 'LUNG700ICRU'; +material{2,2} = -974; +material{2,3} = -724; +material{2,4} = 0.044; +material{2,5} = 0.302; +material{3,1} = 'ICRUTISSUE700ICRU'; +material{3,2} = -724; +material{3,3} = 101; +material{3,4} = 0.302; +material{3,5} = 1.101; +material{4,1} = 'ICRPBONE700ICRU'; +material{4,2} = 101; +material{4,3} = 1976; +material{4,4} = 1.101; +material{4,5} = 2.088; + +% conversion from HU to densities & materials +for s = 1:ct.numOfCtScen + + HUcube{s} = matRad_interp3(dij.ctGrid.x,dij.ctGrid.y',dij.ctGrid.z,ct.cubeHU{s}, ... + dij.doseGrid.x,dij.doseGrid.y',dij.doseGrid.z,'nearest'); + + % projecting out of bounds HU values where necessary + if max(HUcube{s}(:)) > material{end,3} + matRad_cfg.dispWarning('Projecting out of range HU values'); + HUcube{s}(HUcube{s}(:) > material{end,3}) = material{end,3}; + end + if min(HUcube{s}(:)) < material{1,2} + matRad_cfg.dispWarning('Projecting out of range HU values'); + HUcube{s}(HUcube{s}(:) < material{1,2}) = material{1,2}; + end + + % find material index + cubeMatIx{s} = NaN*ones(dij.doseGrid.dimensions,'int32'); + for i = size(material,1):-1:1 + cubeMatIx{s}(HUcube{s} <= material{i,3}) = i; + end + + % create an artificial HU lookup table + hlut = []; + for i = 1:size(material,1) + hlut = [hlut;material{i,2} material{i,4};material{i,3}-1e-10 material{i,5}]; % add eps for interpolation + end + + cubeRho{s} = interp1(hlut(:,1),hlut(:,2),HUcube{s}); + +end + +ompMCgeo.material = material; + +scale = 10; % to convert to cm + +ompMCgeo.xBounds = (dij.doseGrid.resolution.y * (0.5 + [0:dij.doseGrid.dimensions(1)])) ./ scale; +ompMCgeo.yBounds = (dij.doseGrid.resolution.x * (0.5 + [0:dij.doseGrid.dimensions(2)])) ./ scale; +ompMCgeo.zBounds = (dij.doseGrid.resolution.z * (0.5 + [0:dij.doseGrid.dimensions(3)])) ./ scale; + +%% debug visualization +if visBool + + figure + hold on + + axis equal + + % ct box + ctCorner1 = [ompMCgeo.xBounds(1) ompMCgeo.yBounds(1) ompMCgeo.zBounds(1)]; + ctCorner2 = [ompMCgeo.xBounds(end) ompMCgeo.yBounds(end) ompMCgeo.zBounds(end)]; + plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner1(3)],'k' ) + plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) + plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) + plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner1(3) ctCorner1(3)],'k' ) + plot3([ctCorner1(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner2(3) ctCorner2(3)],'k' ) + plot3([ctCorner1(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) + plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) + plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner2(2)],[ctCorner2(3) ctCorner2(3)],'k' ) + plot3([ctCorner1(1) ctCorner1(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' ) + plot3([ctCorner2(1) ctCorner2(1)],[ctCorner1(2) ctCorner1(2)],[ctCorner1(3) ctCorner2(3)],'k' ) + plot3([ctCorner1(1) ctCorner1(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' ) + plot3([ctCorner2(1) ctCorner2(1)],[ctCorner2(2) ctCorner2(2)],[ctCorner1(3) ctCorner2(3)],'k' ) + + xlabel('x [cm]') + ylabel('y [cm]') + zlabel('z [cm]') + + rotate3d on + +end + +%ompMC for matRad returns dose/history * nHistories. +% This factor calibrates to 1 Gy in a %(5x5)cm^2 open field (1 bixel) at +% 5cm depth for SSD = 900 which corresponds to the calibration for the +% analytical base data. +absCalibrationFactor = 3.49056 * 1e12; %Approximate! + +%Now we have to calibrate to the the beamlet width. +absCalibrationFactor = absCalibrationFactor * (pln.propStf.bixelWidth/50)^2; + +outputVariance = matRad_cfg.propMC.ompMC_defaultOutputVariance; + +if isfield(pln,'propMC') && isfield(pln.propMC,'outputVariance') + outputVariance = pln.propMC.outputVariance; +end + +%% Create beamlet source +useCornersSCD = true; %false -> use ISO corners + +scenCount = 0; + +for shiftScen = 1:pln.multScen.totNumShiftScen + + % manipulate isocenter + for k = 1:length(stf) + stf(k).isoCenter = stf(k).isoCenter + pln.multScen.isoShift(shiftScen,:); + end + + for ctScen = 1:pln.multScen.numOfCtScen + for rangeShiftScen = 1:pln.multScen.totNumRangeScen + if pln.multScen.scenMask(ctScen,shiftScen,rangeShiftScen) + scenCount = scenCount + 1; + + numOfBixels = [stf(:).numOfRays]; + beamSource = zeros(dij.numOfBeams, 3); + + bixelCorner = zeros(dij.totalNumOfBixels,3); + bixelSide1 = zeros(dij.totalNumOfBixels,3); + bixelSide2 = zeros(dij.totalNumOfBixels,3); + + counter = 0; + + for i = 1:dij.numOfBeams % loop over all beams + + % define beam source in physical coordinate system in cm + beamSource(i,:) = (stf(i).sourcePoint + stf(i).isoCenter)/10; + + for j = 1:stf(i).numOfRays % loop over all rays / for photons we only have one bixel per ray! + + counter = counter + 1; + + dij.beamNum(counter) = i; + dij.rayNum(counter) = j; + dij.bixelNum(counter) = j; + + if useCornersSCD + beamletCorners = stf(i).ray(j).rayCorners_SCD; + else + beamletCorners = stf(i).ray(j).beamletCornersAtIso; + end + + % get bixel corner and delimiting vectors. + % a) change coordinate system (Isocenter cs-> physical cs) and units mm -> cm + currCorner = (beamletCorners(1,:) + stf(i).isoCenter) ./ scale; + bixelCorner(counter,:) = currCorner; + bixelSide1(counter,:) = (beamletCorners(2,:) + stf(i).isoCenter) ./ scale - currCorner; + bixelSide2(counter,:) = (beamletCorners(4,:) + stf(i).isoCenter) ./ scale - currCorner; + + if visBool + for k = 1:4 + currCornerVis = (beamletCorners(k,:) + stf(i).isoCenter)/10; + % rays connecting source and ray corner + plot3([beamSource(i,1) currCornerVis(1)],[beamSource(i,2) currCornerVis(2)],[beamSource(i,3) currCornerVis(3)],'b') + % connection between corners + lRayCorner = (beamletCorners(mod(k,4) + 1,:) + stf(i).isoCenter)/10; + plot3([lRayCorner(1) currCornerVis(1)],[lRayCorner(2) currCornerVis(2)],[lRayCorner(3) currCornerVis(3)],'r') + end + end + + end + + end + + ompMCsource.nBeams = dij.numOfBeams; + ompMCsource.iBeam = dij.beamNum(:); + + % Switch x and y directions to match ompMC cs. + ompMCsource.xSource = beamSource(:,2); + ompMCsource.ySource = beamSource(:,1); + ompMCsource.zSource = beamSource(:,3); + + ompMCsource.nBixels = sum(numOfBixels(:)); + ompMCsource.xCorner = bixelCorner(:,2); + ompMCsource.yCorner = bixelCorner(:,1); + ompMCsource.zCorner = bixelCorner(:,3); + + ompMCsource.xSide1 = bixelSide1(:,2); + ompMCsource.ySide1 = bixelSide1(:,1); + ompMCsource.zSide1 = bixelSide1(:,3); + + ompMCsource.xSide2 = bixelSide2(:,2); + ompMCsource.ySide2 = bixelSide2(:,1); + ompMCsource.zSide2 = bixelSide2(:,3); + + if visBool + plot3(ompMCsource.ySource,ompMCsource.xSource,ompMCsource.zSource,'rx') + end + + %% Call the OmpMC interface + if pln.multScen.totNumScen > 1 + matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation... \n'); + else + matRad_cfg.dispInfo('matRad: OmpMC photon dose calculation for scenario %d of %d... \n',scenCount,pln.multScen.totNumScen); + end + + ompMCgeo.isoCenter = [stf(:).isoCenter]; + + %Call the Monte Carlo simulation and catch possible mex + %interface issues + try + %If we ask for variance, a field in the dij will be filled + if outputVariance + [dij.physicalDose{ctScen,shiftScen,rangeShiftScen},dij.physicalDose_MCvar{ctScen,shiftScen,rangeShiftScen}] = omc_matrad(cubeRho{ctScen},cubeMatIx{ctScen},ompMCgeo,ompMCsource,ompMCoptions); + else + [dij.physicalDose{ctScen,shiftScen,rangeShiftScen}] = omc_matrad(cubeRho{ctScen},cubeMatIx{ctScen},ompMCgeo,ompMCsource,ompMCoptions); + end + catch ME + errorString = [ME.message '\nThis error was thrown by the MEX-interface of ompMC.\nMex interfaces can raise compatability issues which may be resolved by compiling them by hand directly on your particular system.']; + matRad_cfg.dispError(errorString); + end + + dij.physicalDose{ctScen,shiftScen,rangeShiftScen} = dij.physicalDose{ctScen,shiftScen,rangeShiftScen} * absCalibrationFactor; + if isfield(dij,'physicalDose_MCvar') + dij.physicalDose_MCvar{s} = dij.physicalDose_MCvar{ctScen,shiftScen,rangeShiftScen} * absCalibrationFactor^2; + end + + ompMCgeo.isoCenter = [stf(:).isoCenter]; + + matRad_cfg.dispInfo('matRad: MC photon dose calculation done!\n'); + + try + % wait 0.1s for closing all waitbars + allWaitBarFigures = findall(0,'type','figure','tag','TMWWaitbar'); + delete(allWaitBarFigures); + pause(0.1); + catch + end + + % manipulate isocenter back + for k = 1:length(stf) + stf(k).isoCenter = stf(k).isoCenter - pln.multScen.isoShift(shiftScen,:); + end + end + end + end +end + + +end diff --git a/topas/MatRad_TopasConfig.m b/topas/MatRad_TopasConfig.m index d7ecda0d0..602e0aada 100644 --- a/topas/MatRad_TopasConfig.m +++ b/topas/MatRad_TopasConfig.m @@ -33,7 +33,7 @@ %Simulation parameters numThreads = 0; %number of used threads, 0 = max number of threads (= num cores) - numOfRuns = 5; %Default number of runs / batches + numOfRuns = 1; %Default number of runs / batches modeHistories = 'num'; %'frac'; fracHistories = 1e-4; %Fraction of histories to compute @@ -68,11 +68,8 @@ 'loadConverterFromFile',false); % set true if you want to use your own SchneiderConverter written in "TOPAS_SchneiderConverter" arrayOrdering = 'F'; %'C'; - rsp_basematerial = 'Water'; - rsp_vecLength = 10000; %Bins for RSP values when using RSP conversion with custom image converter - rsp_methodCube = 2; %1: TsBox with variable voxels, 2: TsImageCube with Density Bins and Custom Image converter - + %Scoring scorer = struct('volume',false,... 'doseToMedium',true,... @@ -99,7 +96,7 @@ radiationMode; modules_protons = {'g4em-standard_opt4','g4h-phy_QGSP_BIC_HP','g4decay','g4h-elastic_HP','g4stopping','g4ion-QMD','g4radioactivedecay'}; modules_GenericIon = {'g4em-standard_opt4','g4h-phy_QGSP_BIC_HP','g4decay','g4h-elastic_HP','g4stopping','g4ion-QMD','g4radioactivedecay'}; - modules_gamma = {'g4em-standard_opt4','g4h-phy_QGSP_BIC_HP','g4decay'}; + modules_photons = {'g4em-standard_opt4','g4h-phy_QGSP_BIC_HP','g4decay'}; %Geometry / World worldMaterial = 'G4_AIR'; @@ -133,9 +130,15 @@ 'Scorer_RBE_LEM1','TOPAS_scorer_doseRBE_LEM1.txt.in',... 'Scorer_RBE_WED','TOPAS_scorer_doseRBE_Wedenberg.txt.in',... 'Scorer_RBE_MCN','TOPAS_scorer_doseRBE_McNamara.txt.in'); + + infilenames_photon = struct('virtualGaussian','beamSetup/TOPAS_beamSetup_virtualGaussian.txt.in',... + 'phasespace','beamSetup/TOPAS_beamSetup_Phasespace.txt.in',... + 'uniform','beamSetup/TOPAS_beamSetup_Uniform.txt.in',... + 'mlc','beamSetup/TOPAS_beamSetup_mlc.txt.in' ); + end - properties (SetAccess = private) + properties% (SetAccess = private) thisFolder; MCparam; %Struct with parameters of last simulation to be saved to file @@ -143,18 +146,17 @@ methods function obj = MatRad_TopasConfig() - % MatRad_TopasConfig Construct configuration Class for TOPAS matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class + % Default execution paths are set here + obj.thisFolder = fileparts(mfilename('fullpath')); + obj.workingDir = [obj.thisFolder filesep 'MCrun' filesep]; + % Set default histories from MatRad_Config if isfield(matRad_cfg.propMC,'defaultNumHistories') obj.numHistories = matRad_cfg.propMC.defaultNumHistories; end - % Default execution paths are set here - obj.thisFolder = fileparts(mfilename('fullpath')); - obj.workingDir = [obj.thisFolder filesep 'MCrun' filesep]; - %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) @@ -168,6 +170,405 @@ end end + + function writeStfFieldsPhotons(obj,ct,stf,pln,w) + matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class + + %snippet for future dij simulation + %if ~exist('w','var') + % numBixels = sum([stf(:).totalNumOfBixels)]; + % w = ones(numBixels,1); + %end + + %Bookkeeping + obj.MCparam.nbFields = length(stf); + + %Sanity check + if numel(w) ~= sum([stf(:).totalNumOfBixels]) + matRad_cfg.dispError('Given number of weights (#%d) doesn''t match bixel count in stf (#%d)',numel(w), sum([stf(:).totalNumOfBixels])); + end + + + + nParticlesTotalBixel = round(1e6*w); + nParticlesTotal = sum(nParticlesTotalBixel); + maxParticlesBixel = 1e6*max(w(:)); + minParticlesBixel = round(max([obj.minRelWeight*maxParticlesBixel,1])); + + switch obj.modeHistories + case 'num' + obj.fracHistories = obj.numHistories ./ sum(nParticlesTotalBixel); + case 'frac' + obj.numHistories = sum(nParticlesTotalBixel); + otherwise + matRad_cfg.dispError('Invalid history setting!'); + end + + nParticlesTotal = 0; + + %Preread beam setup + switch obj.beamProfile + case 'phasespace' + fname = fullfile(obj.thisFolder,obj.infilenames_photon.phasespace); + TOPAS_beamSetup = fileread(fname); + obj.pencilBeamScanning = 0 ; + matRad_cfg.dispInfo('Reading ''%s'' Beam Characteristics from ''%s''\n',obj.beamProfile,fname); + case 'virtualGaussian' + fname = fullfile(obj.thisFolder,obj.infilenames_photon.virtualGaussian); + TOPAS_beamSetup = fileread(fname); + matRad_cfg.dispInfo('Reading ''%s'' Beam Characteristics from ''%s''\n',obj.beamProfile,fname) + case 'uniform' + fname = fullfile(obj.thisFolder,obj.infilenames_photon.uniform); + TOPAS_beamSetup = fileread(fname); + matRad_cfg.dispInfo('Reading ''%s'' Beam Characteristics from ''%s''\n',obj.beamProfile,fname); + otherwise + matRad_cfg.dispError('Beam Type ''%s'' not supported for photons',obj.beamProfile); + end + + if isfield(stf(1).ray, 'shapes') + fname = fullfile(obj.thisFolder,obj.infilenames_photon.mlc); + TOPAS_mlcSetup = fileread(fname); + end + + for beamIx = 1:length(stf) + + SAD = stf(beamIx).SAD; + + sourceToNozzleDistance = 0; + nozzleToAxisDistance = SAD; + + %get beamlet properties for each bixel in the stf and write + %it into dataTOPAS + currentBixel = 1; + cutNumOfBixel = 0; + nBeamParticlesTotal(beamIx) = 0; + + collectBixelIdx = []; + + %Loop over rays and then over spots on ray + for rayIx = 1:stf(beamIx).numOfRays + for bixelIx = 1:stf(beamIx).numOfBixelsPerRay(rayIx) + + nCurrentParticles = nParticlesTotalBixel(currentBixel); + + % check whether there are (enough) particles for beam delivery + if (nCurrentParticles>minParticlesBixel) + + collectBixelIdx(end+1) = bixelIx; + cutNumOfBixel = cutNumOfBixel + 1; + bixelEnergy = stf(beamIx).ray(rayIx).energy; + + voxel_x = -stf(beamIx).ray(rayIx).rayPos_bev(3); + voxel_y = stf(beamIx).ray(rayIx).rayPos_bev(1); + + dataTOPAS(cutNumOfBixel).posX = -1.*voxel_x; + dataTOPAS(cutNumOfBixel).posY = voxel_y; + + dataTOPAS(cutNumOfBixel).current = uint32(obj.fracHistories*nCurrentParticles / obj.numOfRuns); + + if obj.pencilBeamScanning + % angleX corresponds to the rotation around the X axis necessary to move the spot in the Y direction + % angleY corresponds to the rotation around the Y' axis necessary to move the spot in the X direction + % note that Y' corresponds to the Y axis after the rotation of angleX around X axis + dataTOPAS(cutNumOfBixel).angleX = atan(dataTOPAS(cutNumOfBixel).posY / SAD); + dataTOPAS(cutNumOfBixel).angleY = atan(-dataTOPAS(cutNumOfBixel).posX ./ (SAD ./ cos(dataTOPAS(cutNumOfBixel).angleX))); + dataTOPAS(cutNumOfBixel).posX = (dataTOPAS(cutNumOfBixel).posX / SAD)*(SAD-nozzleToAxisDistance); + dataTOPAS(cutNumOfBixel).posY = (dataTOPAS(cutNumOfBixel).posY / SAD)*(SAD-nozzleToAxisDistance); + end + + dataTOPAS(cutNumOfBixel).energy = bixelEnergy; + dataTOPAS(cutNumOfBixel).energySpread = 0; + + nBeamParticlesTotal(beamIx) = nBeamParticlesTotal(beamIx) + nCurrentParticles; + + + end + + currentBixel = currentBixel + 1; + + end + end + + + nParticlesTotal = nParticlesTotal + nBeamParticlesTotal(beamIx); + + % discard data if the current has unphysical values + idx = find([dataTOPAS.current] < 1); + dataTOPAS(idx) = []; + cutNumOfBixel = length(dataTOPAS(:)); + + historyCount(beamIx) = uint32(obj.fracHistories * nBeamParticlesTotal(beamIx) / obj.numOfRuns); + + if historyCount(beamIx) < cutNumOfBixel || cutNumOfBixel == 0 + matRad_cfg.dispError('Insufficient number of histories!') + end + + while sum([dataTOPAS.current]) ~= historyCount(beamIx) + % Randomly pick an index with the weigth given by the current + idx = 1:length(dataTOPAS); + % Note: as of Octave version 5.1.0, histcounts is not yet implemented + % using histc instead for compatibility with MATLAB and Octave + %[~,~,R] = histcounts(rand(1),cumsum([0;double(transpose([dataTOPAS(:).current]))./double(sum([dataTOPAS(:).current]))])); + [~,R] = histc(rand(1),cumsum([0;double(transpose([dataTOPAS(:).current]))./double(sum([dataTOPAS(:).current]))])); + randIx = idx(R); + + if (sum([dataTOPAS(:).current]) > historyCount(beamIx)) + if dataTOPAS(randIx).current > 1 + dataTOPAS(randIx).current = dataTOPAS(randIx).current-1; + end + else + dataTOPAS(randIx).current = dataTOPAS(randIx).current+1; + end + end + + historyCount(beamIx) = historyCount(beamIx) * obj.numOfRuns; + + + %sort dataTOPAS according to energy + [~,ixSorted] = sort([dataTOPAS(:).energy]); + dataTOPAS = dataTOPAS(ixSorted); + + %write TOPAS data base file + fieldSetupFileName = sprintf('beamSetup_%s_field%d.txt',obj.label,beamIx); + fileID = fopen(fullfile(obj.workingDir,fieldSetupFileName),'w'); + obj.writeFieldHeader(fileID); + + + %Write modality specific info + fprintf(fileID,'s:Sim/ParticleName = "gamma"\n'); + fprintf(fileID,'u:Sim/ParticleMass = 0\n'); + + particleA = 0; + particleZ = 0; + + obj.modules_photons; + if obj.pencilBeamScanning + fprintf(fileID,'d:Sim/GantryAngle = %f deg\n',stf(beamIx).gantryAngle); + fprintf(fileID,'d:Sim/CouchAngle = %f deg\n',stf(beamIx).couchAngle); + + fprintf(fileID,'d:Tf/TimelineStart = 0. ms\n'); + fprintf(fileID,'d:Tf/TimelineEnd = %i ms\n', 10 * cutNumOfBixel); + fprintf(fileID,'i:Tf/NumberOfSequentialTimes = %i\n', cutNumOfBixel); + fprintf(fileID,'dv:Tf/Beam/Spot/Times = %i ', cutNumOfBixel); + fprintf(fileID,num2str(linspace(10,cutNumOfBixel*10,cutNumOfBixel))); + fprintf(fileID,' ms\n'); + %fprintf(fileID,'uv:Tf/Beam/Spot/Values = %i %s\n',cutNumOfBixel,num2str(collectBixelIdx)); + fprintf(fileID,'s:Tf/Beam/Energy/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/Energy/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/Energy/Values = %i ', cutNumOfBixel); + + fprintf(fileID,num2str([dataTOPAS.energy])); %Transform total energy with atomic number + fprintf(fileID,' MeV\n'); + end + + switch obj.beamProfile + case 'biGaussian' + fprintf(fileID,'s:Tf/Beam/EnergySpread/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/EnergySpread/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/EnergySpread/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.energySpread])); + fprintf(fileID,'\n'); + + fprintf(fileID,'s:Tf/Beam/SigmaX/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaX/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaX/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.spotSize])); + fprintf(fileID,' mm\n'); + fprintf(fileID,'s:Tf/Beam/SigmaXPrime/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaXPrime/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/SigmaXPrime/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.divergence])); + fprintf(fileID,'\n'); + fprintf(fileID,'s:Tf/Beam/CorrelationX/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/CorrelationX/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/CorrelationX/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.correlation])); + fprintf(fileID,'\n'); + + fprintf(fileID,'s:Tf/Beam/SigmaY/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaY/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaY/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.spotSize])); + fprintf(fileID,' mm\n'); + fprintf(fileID,'s:Tf/Beam/SigmaYPrime/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/SigmaYPrime/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/SigmaYPrime/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.divergence])); + fprintf(fileID,'\n'); + fprintf(fileID,'s:Tf/Beam/CorrelationY/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/CorrelationY/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/CorrelationY/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.correlation])); + fprintf(fileID,'\n'); + case 'virtualGaussian' + fprintf(fileID,'s:Tf/Beam/EnergySpread/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/EnergySpread/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/EnergySpread/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.energySpread])); + fprintf(fileID,'\n'); + + if isfield(pln.propStf, 'collimation') + %Use field width for now + fprintf(fileID,'d:So/PencilBeam/BeamPositionSpreadX = %d mm\n', pln.propStf.collimation.fieldWidth); + fprintf(fileID,'d:So/PencilBeam/BeamPositionSpreadY = %d mm\n', pln.propStf.collimation.fieldWidth); + + else + %Set some default value + fprintf(fileID,'d:So/PencilBeam/BeamPositionSpreadX = %d mm\n', 30); + fprintf(fileID,'d:So/PencilBeam/BeamPositionSpreadY = %d mm\n', 30); + end + + case 'uniform' + fprintf(fileID,'s:Tf/Beam/EnergySpread/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/EnergySpread/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'uv:Tf/Beam/EnergySpread/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.energySpread])); + fprintf(fileID,'\n'); + + if isfield(pln.propStf, 'collimation') + %Use field width for now + fprintf(fileID,'d:So/PencilBeam/BeamPositionCutoffX = %d mm\n', pln.propStf.collimation.fieldWidth/2); + fprintf(fileID,'d:So/PencilBeam/BeamPositionCutoffY = %d mm\n', pln.propStf.collimation.fieldWidth/2); + + else + %Set some default value + fprintf(fileID,'d:So/PencilBeam/BeamPositionCutoffX = %d mm\n', 15); + fprintf(fileID,'d:So/PencilBeam/BeamPositionCutoffY = %d mm\n', 15); + end + + case 'simple' + fprintf(fileID,'s:Tf/Beam/FocusFWHM/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/FocusFWHM/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/FocusFWHM/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.focusFWHM])); + fprintf(fileID,' mm\n'); + case 'phasespace' + fprintf(fileID,'d:Sim/GantryAngle = %f deg\n',stf(1).gantryAngle); %just one beam angle for now + fprintf(fileID,'d:Sim/CouchAngle = %f deg\n',stf(1).couchAngle); + + fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', 100 + ct.cubeDim(3)*ct.resolution.z);%Not sure if this is correct,100 is SSD and probably distance from surface to isocenter needs to be added + end + + if obj.pencilBeamScanning + fprintf(fileID,'s:Tf/Beam/AngleX/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/AngleX/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/AngleX/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.angleX])); + fprintf(fileID,' rad\n'); + fprintf(fileID,'s:Tf/Beam/AngleY/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/AngleY/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/AngleY/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.angleY])); + fprintf(fileID,' rad\n'); + + + fprintf(fileID,'s:Tf/Beam/PosX/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/PosX/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/PosX/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.posX])); + fprintf(fileID,' mm\n'); + fprintf(fileID,'s:Tf/Beam/PosY/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/PosY/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'dv:Tf/Beam/PosY/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.posY])); + fprintf(fileID,' mm\n'); + + fprintf(fileID,'s:Tf/Beam/Current/Function = "Step"\n'); + fprintf(fileID,'dv:Tf/Beam/Current/Times = Tf/Beam/Spot/Times ms\n'); + fprintf(fileID,'iv:Tf/Beam/Current/Values = %i ', cutNumOfBixel); + fprintf(fileID,num2str([dataTOPAS.current].*10)); + fprintf(fileID,'\n\n'); + + + % NozzleAxialDistance + fprintf(fileID,'d:Ge/Nozzle/TransZ = -%f mm\n', nozzleToAxisDistance); + fprintf(fileID,'d:Ge/Nozzle/RotX = Tf/Beam/AngleX/Value rad\n'); + fprintf(fileID,'d:Ge/Nozzle/RotY = Tf/Beam/AngleY/Value rad\n'); + fprintf(fileID,'d:Ge/Nozzle/RotZ = 0.0 rad\n'); + + end + + fprintf(fileID,'%s\n',TOPAS_beamSetup); + + if exist('TOPAS_mlcSetup') + fprintf(fileID,'%s\n',TOPAS_mlcSetup); + %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 ', 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 ',numOfLeaves); + for i=1:numOfLeaves + fprintf( fileID,'Tf/LeafXMinus%i/Value ',i); + end + fprintf(fileID,'mm \n'); + + for i=1:numOfLeaves + fprintf(fileID,'s:Tf/LeafXMinus%i/Function = "Step"\n',i); + 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,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,leafTimes); + fprintf(fileID,num2str([1:leafTimes]*10,' % 2d')); + fprintf(fileID,' ms\n'); + fprintf(fileID,'dv:Tf/LeafXPlus%i/Values = %i ', i,leafTimes); + fprintf(fileID,num2str(rightLeafPos(i,:),' % 2d')); + fprintf(fileID,' mm\n\n'); + + end + end + + %translate patient according to beam isocenter + fprintf(fileID,'d:Ge/Patient/TransX = %f mm\n',0.5*ct.resolution.x*(ct.cubeDim(2)+1)-stf(beamIx).isoCenter(1)); + fprintf(fileID,'d:Ge/Patient/TransY = %f mm\n',0.5*ct.resolution.y*(ct.cubeDim(1)+1)-stf(beamIx).isoCenter(2)); + fprintf(fileID,'d:Ge/Patient/TransZ = %f mm\n',0.5*ct.resolution.z*(ct.cubeDim(3)+1)-stf(beamIx).isoCenter(3)); + fprintf(fileID,'d:Ge/Patient/RotX=0. deg\n'); + fprintf(fileID,'d:Ge/Patient/RotY=0. deg\n'); + fprintf(fileID,'d:Ge/Patient/RotZ=0. deg\n'); + + %load topas modules depending on the particle type + fprintf(fileID,'# MODULES\n'); + modules = obj.modules_photons; + moduleString = cellfun(@(s) sprintf('"%s"',s),modules,'UniformOutput',false); + fprintf(fileID,'sv:Ph/Default/Modules = %d %s\n',length(modules),strjoin(moduleString,' ')); + + fclose(fileID); + %write run scripts for TOPAS + for runIx = 1:obj.numOfRuns + runFileName = sprintf('%s_field%d_run%d.txt',obj.label,beamIx,runIx); + fileID = fopen(fullfile(obj.workingDir,runFileName),'w'); + obj.writeRunHeader(fileID,beamIx,runIx); + fprintf(fileID,'includeFile = ./%s\n',fieldSetupFileName); + obj.writeScorers(fileID); + fclose(fileID); + end + + + %Bookkeeping + obj.MCparam.nbParticlesTotal = nParticlesTotal; + obj.MCparam.nbHistoriesTotal = sum(historyCount); + obj.MCparam.nbParticlesField = nBeamParticlesTotal; + obj.MCparam.nbHistoriesField = historyCount; + end + end + + function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % constructor to write all TOPAS fils for local or external simulation % @@ -294,8 +695,13 @@ function writeAllFiles(obj,ct,cst,pln,stf,machine,w) % Generate baseData using the MCemittanceBaseData constructor % Write TOPAS beam properties - topasBaseData = MatRad_MCemittanceBaseData(machine,stf); - obj.writeStfFields(ct,stf,topasBaseData,w); + if ~strcmp(machine.meta.radiationMode,'photons') + topasBaseData = MatRad_MCemittanceBaseData(machine,stf); + obj.writeStfFields(ct,stf,w,topasBaseData); + else + obj.writeStfFieldsPhotons(ct,stf,pln,w); + end + % Save simulation parameters to folder obj.writeMCparam(); @@ -893,10 +1299,17 @@ function writeFieldHeader(obj,fID,ctScen) matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class - fprintf(fID,'u:Sim/HalfValue = %d\n',0.5); - fprintf(fID,'u:Sim/SIGMA2FWHM = %d\n',2.354818); - fprintf(fID,'u:Sim/FWHM2SIGMA = %d\n',0.424661); - fprintf(fID,'\n'); + if ~strcmp(obj.beamProfile,'phasespace') + fprintf(fID,'u:Sim/HalfValue = %d\n',0.5); + fprintf(fID,'u:Sim/SIGMA2FWHM = %d\n',2.354818); + fprintf(fID,'u:Sim/FWHM2SIGMA = %d\n',0.424661); + fprintf(fID,'\n'); + else + fprintf(fID,'u:Sim/HalfValue = %d\n',0.5); + fprintf(fID,'u:Sim/SIGMA2FWHM = %d\n',1.1905); + fprintf(fID,'u:Sim/FWHM2SIGMA = %d\n',0.84); + fprintf(fID,'\n'); + end fprintf(fID,'d:Sim/ElectronProductionCut = %f mm\n',obj.electronProductionCut); fprintf(fID,'s:Sim/WorldMaterial = "%s"\n',obj.worldMaterial); @@ -1126,7 +1539,7 @@ function writeScorers(obj,fID) end end - function writeStfFields(obj,ct,stf,baseData,w) + function writeStfFields(obj,ct,stf,w,baseData) matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class @@ -1390,7 +1803,7 @@ function writeStfFields(obj,ct,stf,baseData,w) particleA = 0; particleZ = 0; - obj.modules_gamma; + obj.modules_photons; %} otherwise matRad_cfg.dispError('Invalid radiation mode %s!',stf.radiationMode) @@ -1630,20 +2043,9 @@ function writePatient(obj,ct,pln) % since at the moment TOPAS does not support ushort % the materials should have indexes between 0 and 32767 % therefore, the maximum length of the vector is 32768 - %answer = '' - %while isempty(answer) - % prompt = 'Length of the imaging vector [2-32768]:'; - % answer = input(prompt) - % if answer > 1 && answer <= 32768 - % vecRSPlength=answer - % else - % answer = '' - % end - %end - + matRad_cfg = MatRad_Config.instance(); %Instance of matRad configuration class medium = obj.rsp_basematerial; - if isequal(obj.arrayOrdering,'C') if matRad_cfg.logLevel > 2 matRad_cfg.dispInfo('Exporting cube in C ordering...\n') diff --git a/topas/beamSetup/TOPAS_beamSetup_Phasespace.txt.in b/topas/beamSetup/TOPAS_beamSetup_Phasespace.txt.in new file mode 100644 index 000000000..447602730 --- /dev/null +++ b/topas/beamSetup/TOPAS_beamSetup_Phasespace.txt.in @@ -0,0 +1,10 @@ +s:So/Phasespace/Type = "PhaseSpace" +s:So/Phasespace/PhaseSpaceFileName = "SIEMENS_PRIMUS_6.0_0.10_15.0x15.0." +s:So/Phasespace/Component = "Nozzle" +#b:So/Example/LimitedAssumeFirstParticleIsNewHistory = "true" +#b:So/Example/LimitedAssumeEveryParticleIsNewHistory = "true" +#b:So/Example/LimitedAssumePhotonIsNewHistory = "true" +i:So/Phasespace/PhaseSpaceMultipleUse = 0 +i:So/Phasespace/NumberOfHistoriesInRun = 5 +i:So/Phasespace/PreCheckShowParticleCountAtInterval = 100000 +b:So/Phasespace/PhaseSpacePreCheck = "False" diff --git a/topas/beamSetup/TOPAS_beamSetup_Uniform.txt.in b/topas/beamSetup/TOPAS_beamSetup_Uniform.txt.in new file mode 100644 index 000000000..9a36fac51 --- /dev/null +++ b/topas/beamSetup/TOPAS_beamSetup_Uniform.txt.in @@ -0,0 +1,15 @@ +s:Ge/BeamSpot/Parent = "Nozzle" +s:Ge/BeamSpot/Type = "Group" +d:Ge/BeamSpot/TransX = Tf/Beam/PosX/Value mm +d:Ge/BeamSpot/TransY = Tf/Beam/PosY/Value mm + +s:So/PencilBeam/BeamParticle = Sim/ParticleName +i:So/PencilBeam/NumberOfHistoriesInRun = Tf/Beam/Current/Value +u:So/PencilBeam/BeamEnergySpread = Tf/Beam/EnergySpread/Value +d:So/PencilBeam/BeamEnergy = Tf/Beam/Energy/Value MeV +s:So/PencilBeam/BeamPositionCutoffShape = "Rectangle" # Rectangle or Ellipse (if Flat or Gaussian) +s:So/PencilBeam/Type = "Beam" # Beam, Isotropic, Emittance or PhaseSpace +s:So/PencilBeam/BeamPositionDistribution = "Flat" # +s:So/PencilBeam/Component = "BeamSpot" +s:So/PencilBeam/BeamAngularDistribution = "None" # None, Flat or Gaussian + diff --git a/topas/beamSetup/TOPAS_beamSetup_biGaussian.txt.in b/topas/beamSetup/TOPAS_beamSetup_biGaussian.txt.in index 9c73912ee..1a7fbdf73 100644 --- a/topas/beamSetup/TOPAS_beamSetup_biGaussian.txt.in +++ b/topas/beamSetup/TOPAS_beamSetup_biGaussian.txt.in @@ -1,3 +1,8 @@ +s:Ge/BeamSpot/Parent = "Nozzle" +s:Ge/BeamSpot/Type = "Group" +d:Ge/BeamSpot/TransX = Tf/Beam/PosX/Value mm +d:Ge/BeamSpot/TransY = Tf/Beam/PosY/Value mm + s:So/PencilBeam/BeamParticle = Sim/ParticleName i:So/PencilBeam/NumberOfHistoriesInRun = Tf/Beam/Current/Value u:So/PencilBeam/BeamEnergySpread = Tf/Beam/EnergySpread/Value diff --git a/topas/beamSetup/TOPAS_beamSetup_generic.txt.in b/topas/beamSetup/TOPAS_beamSetup_generic.txt.in index aa04cbe31..da4458724 100644 --- a/topas/beamSetup/TOPAS_beamSetup_generic.txt.in +++ b/topas/beamSetup/TOPAS_beamSetup_generic.txt.in @@ -1,3 +1,8 @@ +s:Ge/BeamSpot/Parent = "Nozzle" +s:Ge/BeamSpot/Type = "Group" +d:Ge/BeamSpot/TransX = Tf/Beam/PosX/Value mm +d:Ge/BeamSpot/TransY = Tf/Beam/PosY/Value mm + s:So/PencilBeam/BeamParticle = Sim/ParticleName d:So/PencilBeam/BeamEnergy = Tf/Beam/Energy/Value MeV d:So/PencilBeam/BeamPositionSpreadX = Tf/Beam/FocusFWHM/Value mm * Sim/FWHM2SIGMA diff --git a/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in b/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in new file mode 100644 index 000000000..d159d2449 --- /dev/null +++ b/topas/beamSetup/TOPAS_beamSetup_mlc.txt.in @@ -0,0 +1,9 @@ +s:Ge/MultiLeafCollimatorA/Type = "TsMultiLeafCollimator" +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 = 2 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/beamSetup/TOPAS_beamSetup_virtualGaussian.txt.in b/topas/beamSetup/TOPAS_beamSetup_virtualGaussian.txt.in new file mode 100644 index 000000000..f777af0ba --- /dev/null +++ b/topas/beamSetup/TOPAS_beamSetup_virtualGaussian.txt.in @@ -0,0 +1,16 @@ +s:Ge/BeamSpot/Parent = "Nozzle" +s:Ge/BeamSpot/Type = "Group" +d:Ge/BeamSpot/TransX = Tf/Beam/PosX/Value mm +d:Ge/BeamSpot/TransY = Tf/Beam/PosY/Value mm + +s:So/PencilBeam/BeamParticle = Sim/ParticleName +i:So/PencilBeam/NumberOfHistoriesInRun = Tf/Beam/Current/Value +u:So/PencilBeam/BeamEnergySpread = Tf/Beam/EnergySpread/Value +d:So/PencilBeam/BeamEnergy = Tf/Beam/Energy/Value MeV +s:So/PencilBeam/BeamPositionCutoffShape = "Ellipse" # Rectangle or Ellipse (if Flat or Gaussian) +s:So/PencilBeam/Type = "Beam" # Beam, Isotropic, Emittance or PhaseSpace +s:So/PencilBeam/BeamPositionDistribution = "Gaussian" # +s:So/PencilBeam/Component = "BeamSpot" +d:So/PencilBeam/BeamPositionCutoffX = 200 mm +d:So/PencilBeam/BeamPositionCutoffY = 200 mm +s:So/PencilBeam/BeamAngularDistribution = "None" # None, Flat or Gaussian \ No newline at end of file diff --git a/topas/world/TOPAS_matRad_geometry.txt.in b/topas/world/TOPAS_matRad_geometry.txt.in index 6e62ec481..4eb3c522e 100644 --- a/topas/world/TOPAS_matRad_geometry.txt.in +++ b/topas/world/TOPAS_matRad_geometry.txt.in @@ -4,10 +4,8 @@ d:Ge/World/HLY=200.0 cm d:Ge/World/HLZ=200.0 cm b:Ge/World/Invisible = "TRUE" -s:Ge/BeamSpot/Parent = "Nozzle" -s:Ge/BeamSpot/Type = "Group" -d:Ge/BeamSpot/TransX = Tf/Beam/PosX/Value mm -d:Ge/BeamSpot/TransY = Tf/Beam/PosY/Value mm +s:Ge/Nozzle/Parent = "Isocenter" +s:Ge/Nozzle/Type = "Group" s:Ge/Isocenter/Parent = "Gantry" s:Ge/Isocenter/Type = "Group" @@ -28,7 +26,4 @@ d:Ge/Couch/TransY = 0. cm d:Ge/Couch/TransZ = 0. cm d:Ge/Couch/RotX = 0. deg d:Ge/Couch/RotY = -1 * Sim/CouchAngle deg -d:Ge/Couch/RotZ = 0. deg - -s:Ge/Nozzle/Parent = "Isocenter" -s:Ge/Nozzle/Type = "Group" \ No newline at end of file +d:Ge/Couch/RotZ = 0. deg \ No newline at end of file