11classdef matRad_TopasMCEngine < DoseEngines .matRad_MonteCarloEngineAbstract
2- % matRad_TopasMCEngine
2+ % matRad_TopasMCEngine
33 % Implementation of the TOPAS interface for Monte Carlo dose
44 % calculation
55 %
1818 %
1919 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2020
21- properties (Constant )
21+ properties (Constant )
2222 possibleRadiationModes = {' photons' ,' protons' ,' helium' ,' carbon' };
2323 name = ' TOPAS' ;
2424 shortName = ' TOPAS' ;
3838 topasExecCommand ; % Defaults will be set during construction according to TOPAS installation instructions and used system
3939
4040 parallelRuns = false ; % Starts runs in parallel
41-
41+
4242 externalCalculation = ' off' ; % Generates folder for external TOPAS calculation (e.g. on a server)
4343
4444 workingDir ; % working directory for the simulation
157157 ' Scorer_RBE_MCN' ,' TOPAS_scorer_doseRBE_McNamara.txt.in' , ...
158158 ... %PhaseSpace Source
159159 ' phaseSpaceSourcePhotons' ,' VarianClinaciX_6MV_20x20_aboveMLC_w2' );
160-
160+
161161
162162 end
163163
182182 function setDefaults(this )
183183 this .setDefaults @DoseEngines .matRad_MonteCarloEngineAbstract();
184184 matRad_cfg = MatRad_Config .instance(); % Instance of matRad configuration class
185-
185+
186186 this.useGivenEqDensityCube = matRad_cfg .defaults .propDoseCalc .useGivenEqDensityCube ;
187-
187+
188188 % Default execution paths are set here
189189 this.topasFolder = [matRad_cfg .matRadSrcRoot filesep ' doseCalc' filesep ' topas' filesep ];
190190 this.workingDir = [matRad_cfg .primaryUserFolder filesep ' TOPAS' filesep ];
@@ -193,7 +193,7 @@ function setDefaults(this)
193193 mkdir(this .workingDir );
194194 matRad_cfg .dispInfo(' Created TOPAS working directory in userfolder %s\n ' ,this .workingDir );
195195 end
196-
196+
197197
198198 % Let's set some default commands taken from topas installation
199199 % instructions for mac & debain/ubuntu
@@ -465,7 +465,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
465465 return ;
466466 else
467467 end
468-
468+
469469
470470 %% Initialize dose grid and dij
471471
@@ -496,12 +496,12 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
496496 end
497497 end
498498 end
499-
499+
500500 for i = 1 : numel(stf )
501501 if strcmp(stf(i ).radiationMode,' photons' )
502502 stf(i ).ray.energy = stf(i ).ray.energy.* ones(size(w ));
503503 end
504- end
504+ end
505505
506506 % Get photon parameters for RBExDose calculation
507507 if this .calcBioDose
@@ -645,7 +645,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
645645 dij.beamNum = 1 ;
646646 dij.bixelNum = 1 ;
647647 dij.ctGrid = this .ctR .ctGrid ;
648- dij.doseGrid = this .doseGrid ;
648+ dij.doseGrid = this .doseGrid ;
649649 dij.numOfBeams = 1 ;
650650 dij.numOfRaysPerBeam = 1 ;
651651 dij.numOfScenarios = this .multScen .totNumScen ;
@@ -671,7 +671,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
671671
672672 end
673673
674-
674+
675675 function dij = initDoseCalc(this ,ct ,cst ,stf )
676676 dij = this .initDoseCalc @DoseEngines .matRad_MonteCarloEngineAbstract(ct ,cst ,stf );
677677 matRad_cfg = MatRad_Config .instance();
@@ -702,35 +702,35 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
702702 % Allpcate resampled cubes
703703 cubeHUresampled = cell(1 ,ct .numOfCtScen );
704704 cubeResampled = cell(1 ,ct .numOfCtScen );
705-
705+
706706 % Perform resampling to dose grid
707707 for s = 1 : ct .numOfCtScen
708708 cubeHUresampled{s } = matRad_interp3(dij .ctGrid .x , dij .ctGrid .y ' , dij .ctGrid .z ,ct.cubeHU{s }, ...
709709 dij .doseGrid .x ,dij .doseGrid .y ' ,dij .doseGrid .z ,' linear' );
710710 cubeResampled{s } = matRad_interp3(dij .ctGrid .x , dij .ctGrid .y ' , dij .ctGrid .z ,ct.cube{s }, ...
711711 dij .doseGrid .x ,dij .doseGrid .y ' ,dij .doseGrid .z ,' linear' );
712712 end
713-
713+
714714 % Allocate temporary resampled CT
715715 this.ctR = ct ;
716716 this.ctR.cube = cell(1 );
717717 this.ctR.cubeHU = cell(1 );
718-
718+
719719 % Set CT resolution to doseGrid resolution
720720 this.ctR.resolution = dij .doseGrid .resolution ;
721721 this.ctR.cubeDim = dij .doseGrid .dimensions ;
722722 this.ctR.x = dij .doseGrid .x ;
723723 this.ctR.y = dij .doseGrid .y ;
724724 this.ctR.z = dij .doseGrid .z ;
725-
725+
726726 % Write resampled cubes
727727 this.ctR.cubeHU = cubeHUresampled ;
728728 this.ctR.cube = cubeResampled ;
729-
729+
730730 % Set flag for complete resampling
731731 this.ctR.resampled = 1 ;
732732 this.ctR.ctGrid = dij .doseGrid ;
733-
733+
734734 % Save original grid
735735 this.ctR.originalGrid = dij .ctGrid ;
736736 matRad_cfg .dispInfo(' done!\n ' );
@@ -905,7 +905,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
905905 end
906906
907907 % Tally per field
908- if isfield(topasSum ,' Sum' )
908+ if isfield(topasSum ,' Sum' )
909909 if obj .calc4DInterplay || obj .MCparam .numOfCtScen > 1
910910 if strcmp(tname , ' physicalDose' )
911911 topasCube.([' phaseDose_beam' num2str(f )]){ctScen } = topasSum .Sum ;
@@ -1071,7 +1071,7 @@ function writeAllFiles(obj,ct,cst,stf,machine,w)
10711071 % Allocate possible scored quantities
10721072 processedQuantities = {' ' ,' _std' ,' _batchStd' };
10731073 topasCubesTallies = unique(erase(topasCubesTallies ,processedQuantities(2 : end )));
1074-
1074+
10751075
10761076 % Loop through 4D scenarios
10771077 for ctScen = 1 : dij .numOfScenarios
@@ -1232,7 +1232,7 @@ function writeRunHeader(obj,fID,fieldIx,runIx,ctScen)
12321232 fprintf(fID ,' \n ' );
12331233 fprintf(fID ,[' i:Ts/Seed = ' ,num2str(runIx ),' \n ' ]);
12341234
1235- % TODO: remove or document
1235+ % TODO: remove or document
12361236 % fprintf(fID,'includeFile = %s/TOPAS_Simulation_Setup.txt\n',obj.thisFolder);
12371237 % fprintf(fID,'includeFile = %s/TOPAS_matRad_geometry.txt\n',obj.thisFolder);
12381238 % fprintf(fID,'includeFile = %s/TOPAS_scorer_surfaceIC.txt\n',obj.thisFolder);
@@ -1295,12 +1295,12 @@ function writeScorers(obj,fID,beamIx)
12951295 matRad_cfg .dispDebug(' Reading doseToMedium scorer from %s\n ' ,fname );
12961296 scorerName = fileread(fname );
12971297 fprintf(fID ,' \n%s\n\n ' ,scorerName );
1298-
1298+
12991299 if obj .calc4DInterplay
13001300 for PhaseNum = obj.MCparam.Phases{beamIx }'
13011301 scorerTxt = fileread(fname );
1302- scorerTxt = replace (scorerTxt , ' /Patient' , [' /Patient' num2str(PhaseNum )]);
1303- scorerTxt = replace (scorerTxt , ' _physicalDose' , [' _physicalDose-matRad_cube' num2str(PhaseNum )]);
1302+ scorerTxt = strrep (scorerTxt , ' /Patient' , [' /Patient' num2str(PhaseNum )]);
1303+ scorerTxt = strrep (scorerTxt , ' _physicalDose' , [' _physicalDose-matRad_cube' num2str(PhaseNum )]);
13041304 fprintf(fID ,' \n%s\n\n ' ,scorerTxt );
13051305 fprintf(fID , ' b:Sc/Patient%i /Tally_DoseToMedium/Active = Tf/Patient%i /Tally_DoseToMedium/Active/Value\n ' , PhaseNum , PhaseNum );
13061306 fprintf(fID ,' s:Tf/Patient%i /Tally_DoseToMedium/Active/Function = "Step"\n ' , PhaseNum );
@@ -1354,33 +1354,33 @@ function writeScorers(obj,fID,beamIx)
13541354 fprintf(fID ,' \n%s\n\n ' ,scorerName );
13551355
13561356 if obj .calc4DInterplay
1357- for PhaseNum = obj.MCparam.Phases{beamIx }'
1357+ for PhaseNum = obj.MCparam.Phases{beamIx }'
13581358 scorerTxt = fileread(fname );
1359- scorerTxt = replace (scorerTxt , ' Alpha/' , [' Alpha' num2str(PhaseNum ) ' /' ]);
1360- scorerTxt = replace (scorerTxt , ' Beta/' , [' Beta' num2str(PhaseNum ) ' /' ]);
1361- scorerTxt = replace (scorerTxt , ' /RBE' , [' /RBE' num2str(PhaseNum )]);
1359+ scorerTxt = strrep (scorerTxt , ' Alpha/' , [' Alpha' num2str(PhaseNum ) ' /' ]);
1360+ scorerTxt = strrep (scorerTxt , ' Beta/' , [' Beta' num2str(PhaseNum ) ' /' ]);
1361+ scorerTxt = strrep (scorerTxt , ' /RBE' , [' /RBE' num2str(PhaseNum )]);
13621362 if contains(lower(obj.scorer.RBE_model{i }),' mcn' )
1363- scorerTxt = replace (scorerTxt , ' _alpha_MCN' , [' _alpha_MCN-matRad_cube' num2str(PhaseNum )]);
1364- scorerTxt = replace (scorerTxt , ' _beta_MCN' , [' _beta_MCN-matRad_cube' num2str(PhaseNum )]);
1363+ scorerTxt = strrep (scorerTxt , ' _alpha_MCN' , [' _alpha_MCN-matRad_cube' num2str(PhaseNum )]);
1364+ scorerTxt = strrep (scorerTxt , ' _beta_MCN' , [' _beta_MCN-matRad_cube' num2str(PhaseNum )]);
13651365 elseif contains(lower(obj.scorer.RBE_model{i }),' wed' )
1366- scorerTxt = replace (scorerTxt , ' _alpha_WED' , [' _alpha_WED-matRad_cube' num2str(PhaseNum )]);
1367- scorerTxt = replace (scorerTxt , ' _beta_WED' , [' _beta_WED-matRad_cube' num2str(PhaseNum )]);
1366+ scorerTxt = strrep (scorerTxt , ' _alpha_WED' , [' _alpha_WED-matRad_cube' num2str(PhaseNum )]);
1367+ scorerTxt = strrep (scorerTxt , ' _beta_WED' , [' _beta_WED-matRad_cube' num2str(PhaseNum )]);
13681368 elseif contains(lower(obj.scorer.RBE_model{i }),' lem' )
1369- scorerTxt = replace (scorerTxt , ' _alpha_LEM' , [' _alpha_LEM-matRad_cube' num2str(PhaseNum )]);
1370- scorerTxt = replace (scorerTxt , ' _beta_LEM' , [' _beta_LEM-matRad_cube' num2str(PhaseNum )]);
1371- scorerTxt = replace (scorerTxt , ' _RBE_LEM' , [' _RBE_LEM-matRad_cube' num2str(PhaseNum )]);
1369+ scorerTxt = strrep (scorerTxt , ' _alpha_LEM' , [' _alpha_LEM-matRad_cube' num2str(PhaseNum )]);
1370+ scorerTxt = strrep (scorerTxt , ' _beta_LEM' , [' _beta_LEM-matRad_cube' num2str(PhaseNum )]);
1371+ scorerTxt = strrep (scorerTxt , ' _RBE_LEM' , [' _RBE_LEM-matRad_cube' num2str(PhaseNum )]);
13721372 elseif contains(lower(obj.scorer.RBE_model{i }),' libamtrack' )
1373- scorerTxt = replace (scorerTxt , ' _alpha_libamtrack' , [' _alpha_libamtrack-matRad_cube' num2str(PhaseNum )]);
1374- scorerTxt = replace (scorerTxt , ' _beta_libamtrack' , [' _beta_libamtrack-matRad_cube' num2str(PhaseNum )]);
1375- scorerTxt = replace (scorerTxt , ' _RBE_libamtrack' , [' _RBE_libamtrack-matRad_cube' num2str(PhaseNum )]);
1373+ scorerTxt = strrep (scorerTxt , ' _alpha_libamtrack' , [' _alpha_libamtrack-matRad_cube' num2str(PhaseNum )]);
1374+ scorerTxt = strrep (scorerTxt , ' _beta_libamtrack' , [' _beta_libamtrack-matRad_cube' num2str(PhaseNum )]);
1375+ scorerTxt = strrep (scorerTxt , ' _RBE_libamtrack' , [' _RBE_libamtrack-matRad_cube' num2str(PhaseNum )]);
13761376 end
13771377 % dont allways write cell lines
13781378 if contains(scorerTxt ,' ### HCP Tabulated ###' )
13791379 scorerTxt = extractBefore(scorerTxt , ' ### HCP Tabulated ###' );
13801380 end
13811381 fprintf(fID ,' \n%s\n\n ' ,scorerTxt );
13821382 scorerInString = {' tabulatedAlpha' , ' tabulatedBeta' , ' RBE' , ' McNamaraAlpha' , ' McNamaraBeta' , ' WedenbergAlpha' , ' WedenbergBeta' };
1383- for ixWrite = ixToWrite4DInterplay
1383+ for ixWrite = ixToWrite4DInterplay
13841384 fprintf(fID ,' b:Sc/%s%i /Active = Tf/%s%i /Active/Value\n ' , scorerInString{ixWrite },PhaseNum ,scorerInString{ixWrite },PhaseNum );
13851385 fprintf(fID ,' s:Tf/%s%i /Active/Function = "Step"\n ' , scorerInString{ixWrite }, PhaseNum );
13861386 fprintf(fID ,' dv:Tf/%s%i /Active/Times= %i ' ,scorerInString{ixWrite },PhaseNum , obj .MCparam .cutNumOfBixel(beamIx ));
@@ -1450,7 +1450,7 @@ function writeScorers(obj,fID,beamIx)
14501450 end
14511451 fprintf(fID ,' s:Sc/%s%s /ReferencedSubScorer_Dose = "Tally_DoseToWater"\n ' ,scorerPrefix ,scorerNames{s });
14521452 if obj .calc4DInterplay
1453- for PhaseNum = obj.MCparam.Phases{beamIx }'
1453+ for PhaseNum = obj.MCparam.Phases{beamIx }'
14541454 if strcmp(obj .radiationMode ,' protons' )
14551455 fprintf(fID ,' s:Sc/%s%s%i /ReferencedSubScorer_LET = "ProtonLET%i "\n ' ,scorerPrefix ,scorerNames{s },PhaseNum ,PhaseNum );
14561456 end
@@ -1473,8 +1473,8 @@ function writeScorers(obj,fID,beamIx)
14731473 if obj .calc4DInterplay
14741474 for PhaseNum = obj.MCparam.Phases{beamIx }'
14751475 scorerTxt = fileread(fname );
1476- scorerTxt = replace (scorerTxt , ' Tally_DoseToWater' , [' Tally_DoseToWater' num2str(PhaseNum )]);
1477- scorerTxt = replace (scorerTxt , ' _doseToWater' , [' _doseToWater-matRad_cube' num2str(PhaseNum )]);
1476+ scorerTxt = strrep (scorerTxt , ' Tally_DoseToWater' , [' Tally_DoseToWater' num2str(PhaseNum )]);
1477+ scorerTxt = strrep (scorerTxt , ' _doseToWater' , [' _doseToWater-matRad_cube' num2str(PhaseNum )]);
14781478 fprintf(fID ,' \n%s\n\n ' ,scorerTxt );
14791479 fprintf(fID , ' b:Sc/Tally_DoseToWater%i /Active = Tf/Tally_DoseToWater%i /Active/Value\n ' , PhaseNum , PhaseNum );
14801480 fprintf(fID ,' s:Tf/Tally_DoseToWater%i /Active/Function = "Step"\n ' , PhaseNum );
@@ -1502,8 +1502,8 @@ function writeScorers(obj,fID,beamIx)
15021502 if obj .calc4DInterplay
15031503 for PhaseNum = obj.MCparam.Phases{beamIx }'
15041504 scorerTxt = fileread(fname );
1505- scorerTxt = replace (scorerTxt , ' /ProtonLET' , [' /ProtonLET' num2str(PhaseNum )]);
1506- scorerTxt = replace (scorerTxt , ' _LET' , [' _LET-matRad_cube' num2str(PhaseNum )]);
1505+ scorerTxt = strrep (scorerTxt , ' /ProtonLET' , [' /ProtonLET' num2str(PhaseNum )]);
1506+ scorerTxt = strrep (scorerTxt , ' _LET' , [' _LET-matRad_cube' num2str(PhaseNum )]);
15071507 fprintf(fID ,' \n%s\n\n ' ,scorerTxt );
15081508 fprintf(fID , ' b:Sc/ProtonLET%i /Active = Tf/ProtonLET%i /Active/Value\n ' , PhaseNum , PhaseNum );
15091509 fprintf(fID ,' s:Tf/ProtonLET%i /Active/Function = "Step"\n ' , PhaseNum );
@@ -1700,8 +1700,8 @@ function writeStfFields(obj,ct,stf,w,baseData)
17001700 selectedData = [];
17011701 focusIndex = baseData .selectedFocus(baseData .energyIndex );
17021702
1703- scalarFields = {' NominalEnergy' ,' EnergySpread' ,' MeanEnergy' };
1704-
1703+ scalarFields = {' NominalEnergy' ,' EnergySpread' ,' MeanEnergy' };
1704+
17051705 for i = 1 : numel(focusIndex )
17061706 for field = scalarFields
17071707 baseData .monteCarloData(i ).(field{1 }) = ones(1 ,max(focusIndex ))*baseData .monteCarloData(i ).(field{1 });
@@ -1821,7 +1821,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
18211821 nBeamParticlesTotal(beamIx ) = nBeamParticlesTotal(beamIx ) + nCurrentParticles ;
18221822 currentBixel = currentBixel + 1 ;
18231823
1824-
1824+
18251825 end
18261826 end
18271827
@@ -1890,7 +1890,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
18901890 end
18911891
18921892 % NozzleAxialDistance
1893- if isPhoton
1893+ if isPhoton
18941894 fprintf(fileID ,' d:Ge/Nozzle/TransZ = -%f mm\n ' , stf(beamIx ).SCD+ 40 ); % Phasespace hardcorded infront of MLC at SSD 46 cm
18951895 else
18961896 fprintf(fileID ,' d:Ge/Nozzle/TransZ = -%f mm\n ' , nozzleToAxisDistance );
@@ -2082,7 +2082,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
20822082
20832083 fprintf(fileID ,' d:Sim/GantryAngle = %f deg\n ' ,stf(beamIx ).gantryAngle); % just one beam angle for now
20842084 fprintf(fileID ,' d:Sim/CouchAngle = %f deg\n ' ,stf(beamIx ).couchAngle);
2085- % Here the phasespace file is loaded and referenced in the beamSetup file
2085+ % Here the phasespace file is loaded and referenced in the beamSetup file
20862086 if strcmp(obj .externalCalculation ,' write' )
20872087 matRad_cfg .dispWarning([' External calculation and phaseSpace selected, manually place ' obj .infilenames .phaseSpaceSourcePhotons ' .header and ' obj .infilenames .phaseSpaceSourcePhotons ' .phsp into your simulation directory.' ]);
20882088 else
@@ -2091,7 +2091,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
20912091 end
20922092 end
20932093 % phasespaceStr = ['..' filesep 'beamSetup' filesep 'phasespace' filesep phaseSpaceFileName];
2094- % &phasespaceStr = replace (phasespaceStr, '\', '/');
2094+ % &phasespaceStr = strrep (phasespaceStr, '\', '/');
20952095 fprintf(fileID ,' s:So/Phasespace/PhaseSpaceFileName = "%s "\n ' , obj .infilenames .phaseSpaceSourcePhotons );
20962096
20972097 end
@@ -2238,7 +2238,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
22382238 if isempty(obj .calcTimeSequence )
22392239 matRad_cfg .dispError(' Time Sequence Data required for 4D calculation \n ' );
22402240 end
2241-
2241+
22422242 % Write changing CT phases in matRad_cube file
22432243 outfilePatient = obj .outfilenames .patientParam ;
22442244 outfilePatient = strsplit(outfilePatient ,' .' );
@@ -2261,7 +2261,7 @@ function writeStfFields(obj,ct,stf,w,baseData)
22612261 fprintf(fIDPatient ,' sv:Tf/InputFile/Values = %d %s \n ' , numel(fileNamesPhases ),strjoin(fileNamesPhases ,' ' ));
22622262 fprintf(fIDPatient , ' \n ' );
22632263 fclose(fIDPatient );
2264-
2264+
22652265 % write On of Feature into MCparam to create scorers
22662266 % later
22672267 uniquePhaseNum = unique(PhaseNum );
@@ -2714,7 +2714,7 @@ function writeMCparam(obj)
27142714 % Used to check against a machine file if a specific quantity can be
27152715 % computed.
27162716 function q = providedQuantities(machine )
2717- q = {' physicalDose' ,' LET' ,' alpha' ,' beta' };
2717+ q = {' physicalDose' ,' LET' ,' alpha' ,' beta' };
27182718 end
27192719 end
27202720end
0 commit comments